"Fossies" - the Fresh Open Source Software Archive

Member "meep-1.10.0/tests/two_dimensional.cpp" (27 Feb 2019, 13387 Bytes) of package /linux/privat/meep-1.10.0.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. See also the last Fossies "Diffs" side-by-side code changes report for "two_dimensional.cpp": 1.8.0_vs_1.9.0.

    1 /* Copyright (C) 2005-2019 Massachusetts Institute of Technology
    2 %
    3 %  This program is free software; you can redistribute it and/or modify
    4 %  it under the terms of the GNU General Public License as published by
    5 %  the Free Software Foundation; either version 2, or (at your option)
    6 %  any later version.
    7 %
    8 %  This program is distributed in the hope that it will be useful,
    9 %  but WITHOUT ANY WARRANTY; without even the implied warranty of
   10 %  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   11 %  GNU General Public License for more details.
   12 %
   13 %  You should have received a copy of the GNU General Public License
   14 %  along with this program; if not, write to the Free Software Foundation,
   15 %  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
   16 */
   17 
   18 #include <stdio.h>
   19 #include <stdlib.h>
   20 #include <signal.h>
   21 
   22 #include <meep.hpp>
   23 using namespace meep;
   24 using std::complex;
   25 
   26 double one(const vec &) { return 1.0; }
   27 double targets(const vec &pt) {
   28   const double r = sqrt(pt.x() * pt.x() + pt.y() * pt.y());
   29   double dr = r;
   30   while (dr > 1)
   31     dr -= 1;
   32   if (dr > 0.7001) return 12.0;
   33   return 1.0;
   34 }
   35 
   36 #if MEEP_SINGLE
   37 static const double tol = 1e-3, thresh = 1e-5;
   38 #else
   39 static const double tol = 1e-11, thresh = 1e-5;
   40 #endif
   41 
   42 int compare(double a, double b, const char *n) {
   43   if (fabs(a - b) > fabs(b) * tol && fabs(b) > thresh) {
   44     master_printf("%s differs by\t%g out of\t%g\n", n, a - b, b);
   45     master_printf("This gives a fractional error of %g\n", fabs(a - b) / fabs(b));
   46     return 0;
   47   } else {
   48     return 1;
   49   }
   50 }
   51 
   52 int compare_point(fields &f1, fields &f2, const vec &p) {
   53   monitor_point m1, m_test;
   54   f1.get_point(&m_test, p);
   55   f2.get_point(&m1, p);
   56   for (int i = 0; i < 10; i++) {
   57     component c = (component)i;
   58     if (f1.gv.has_field(c)) {
   59       complex<double> v1 = m_test.get_component(c), v2 = m1.get_component(c);
   60       if (abs(v1 - v2) > tol * abs(v2) && abs(v2) > thresh) {
   61         master_printf("%s differs:  %g %g out of %g %g\n", component_name(c), real(v2 - v1),
   62                       imag(v2 - v1), real(v2), imag(v2));
   63         master_printf("This comes out to a fractional error of %g\n", abs(v1 - v2) / abs(v2));
   64         master_printf("Right now I'm looking at %g %g, time %g\n", p.x(), p.y(), f1.time());
   65         return 0;
   66       }
   67     }
   68   }
   69   return 1;
   70 }
   71 
   72 int test_metal(double eps(const vec &), int splitting, const char *mydirname) {
   73   double a = 10.0;
   74   double ttot = 17.0;
   75 
   76   grid_volume gv = voltwo(3.0, 2.0, a);
   77   structure s1(gv, eps);
   78   structure s(gv, eps, no_pml(), identity(), splitting);
   79   s.set_output_directory(mydirname);
   80   s1.set_output_directory(mydirname);
   81 
   82   s.add_susceptibility(one, E_stuff, lorentzian_susceptibility(0.3, 0.1));
   83   s1.add_susceptibility(one, E_stuff, lorentzian_susceptibility(0.3, 0.1));
   84 
   85   master_printf("Metal+dispersion test using %d chunks...\n", splitting);
   86   fields f(&s);
   87   f.add_point_source(Hz, 0.7, 2.5, 0.0, 4.0, vec(0.3, 0.5), 1.0);
   88   f.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299, 0.401), 1.0);
   89   fields f1(&s1);
   90   f1.add_point_source(Hz, 0.7, 2.5, 0.0, 4.0, vec(0.3, 0.5), 1.0);
   91   f1.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299, 0.401), 1.0);
   92   double field_energy_check_time = 8.0;
   93   while (f.time() < ttot) {
   94     f.step();
   95     f1.step();
   96     if (!compare_point(f, f1, vec(0.5, 0.01))) return 0;
   97     if (!compare_point(f, f1, vec(0.46, 0.33))) return 0;
   98     if (!compare_point(f, f1, vec(1.0, 1.0))) return 0;
   99     if (f.time() >= field_energy_check_time) {
  100       if (!compare(f.field_energy(), f1.field_energy(), "   total energy")) return 0;
  101       if (!compare(f.electric_energy_in_box(gv.surroundings()),
  102                    f1.electric_energy_in_box(gv.surroundings()), "electric energy"))
  103         return 0;
  104       if (!compare(f.magnetic_energy_in_box(gv.surroundings()),
  105                    f1.magnetic_energy_in_box(gv.surroundings()), "magnetic energy"))
  106         return 0;
  107       field_energy_check_time += 5.0;
  108     }
  109   }
  110   return 1;
  111 }
  112 
  113 int test_periodic(double eps(const vec &), int splitting, const char *mydirname) {
  114   double a = 10.0;
  115   double ttot = 17.0;
  116 
  117   grid_volume gv = voltwo(3.0, 2.0, a);
  118   structure s1(gv, eps);
  119   structure s(gv, eps, no_pml(), identity(), splitting);
  120   s.set_output_directory(mydirname);
  121   s1.set_output_directory(mydirname);
  122 
  123   master_printf("Periodic test using %d chunks...\n", splitting);
  124   fields f(&s);
  125   f.use_bloch(vec(0.1, 0.7));
  126   f.add_point_source(Hz, 0.7, 2.5, 0.0, 4.0, vec(0.3, 0.5), 1.0);
  127   f.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299, 0.401), 1.0);
  128   fields f1(&s1);
  129   f1.use_bloch(vec(0.1, 0.7));
  130   f1.add_point_source(Hz, 0.7, 2.5, 0.0, 4.0, vec(0.3, 0.5), 1.0);
  131   f1.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299, 0.401), 1.0);
  132   double field_energy_check_time = 8.0;
  133   while (f.time() < ttot) {
  134     f.step();
  135     f1.step();
  136     if (!compare_point(f, f1, vec(0.5, 0.01))) return 0;
  137     if (!compare_point(f, f1, vec(0.46, 0.33))) return 0;
  138     if (!compare_point(f, f1, vec(1.0, 1.0))) return 0;
  139     if (f.time() >= field_energy_check_time) {
  140       if (!compare(f.field_energy(), f1.field_energy(), "   total energy")) return 0;
  141       if (!compare(f.electric_energy_in_box(gv.surroundings()),
  142                    f1.electric_energy_in_box(gv.surroundings()), "electric energy"))
  143         return 0;
  144       if (!compare(f.magnetic_energy_in_box(gv.surroundings()),
  145                    f1.magnetic_energy_in_box(gv.surroundings()), "magnetic energy"))
  146         return 0;
  147       field_energy_check_time += 5.0;
  148     }
  149   }
  150   return 1;
  151 }
  152 
  153 int test_periodic_tm(double eps(const vec &), int splitting, const char *mydirname) {
  154   double a = 10.0;
  155   double ttot = 17.0;
  156 
  157   grid_volume gv = voltwo(3.0, 2.0, a);
  158   structure s1(gv, eps);
  159   structure s(gv, eps, no_pml(), identity(), splitting);
  160   s.set_output_directory(mydirname);
  161   s1.set_output_directory(mydirname);
  162 
  163   master_printf("Periodic 2D TM test using %d chunks...\n", splitting);
  164   fields f(&s);
  165   f.use_bloch(vec(0.1, 0.7));
  166   f.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299, 0.401), 1.0);
  167   fields f1(&s1);
  168   f1.use_bloch(vec(0.1, 0.7));
  169   f1.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(1.299, 0.401), 1.0);
  170   double field_energy_check_time = 8.0;
  171   while (f.time() < ttot) {
  172     f.step();
  173     f1.step();
  174     if (!compare_point(f, f1, vec(0.5, 0.01))) return 0;
  175     if (!compare_point(f, f1, vec(0.46, 0.33))) return 0;
  176     if (!compare_point(f, f1, vec(1.0, 1.0))) return 0;
  177     if (f.time() >= field_energy_check_time) {
  178       if (!compare(f.field_energy(), f1.field_energy(), "   total energy")) return 0;
  179       if (!compare(f.electric_energy_in_box(gv.surroundings()),
  180                    f1.electric_energy_in_box(gv.surroundings()), "electric energy"))
  181         return 0;
  182       if (!compare(f.magnetic_energy_in_box(gv.surroundings()),
  183                    f1.magnetic_energy_in_box(gv.surroundings()), "magnetic energy"))
  184         return 0;
  185       field_energy_check_time += 5.0;
  186     }
  187   }
  188   return 1;
  189 }
  190 
  191 int test_pml(double eps(const vec &), int splitting, const char *mydirname) {
  192   double a = 10.0;
  193 
  194   grid_volume gv = voltwo(3.0, 2.0, a);
  195   structure s1(gv, eps, pml(1.0, X) + pml(1.0, Y, High));
  196   structure s(gv, eps, pml(1.0, X) + pml(1.0, Y, High), identity(), splitting);
  197   s.set_output_directory(mydirname);
  198   s1.set_output_directory(mydirname);
  199 
  200   master_printf("Testing pml while splitting into %d chunks...\n", splitting);
  201   fields f(&s);
  202   f.add_point_source(Hz, 0.7, 1.5, 0.0, 4.0, vec(1.5, 0.5), 1.0);
  203   f.add_point_source(Ez, 0.8, 1.6, 0.0, 4.0, vec(1.299, 0.401), 1.0);
  204   fields f1(&s1);
  205   f1.add_point_source(Hz, 0.7, 1.5, 0.0, 4.0, vec(1.5, 0.5), 1.0);
  206   f1.add_point_source(Ez, 0.8, 1.6, 0.0, 4.0, vec(1.299, 0.401), 1.0);
  207   const double deltaT = 100.0;
  208   const double ttot = 3.1 * deltaT;
  209   double field_energy_check_time = deltaT;
  210 
  211   while (f.time() < f.last_source_time())
  212     f.step();
  213   while (f1.time() < f1.last_source_time())
  214     f1.step();
  215 
  216   double last_energy = f.field_energy();
  217   while (f.time() < ttot) {
  218     f.step();
  219     f1.step();
  220     if (f.time() >= field_energy_check_time) {
  221       if (!compare_point(f, f1, vec(0.5, 0.01))) return 0;
  222       if (!compare_point(f, f1, vec(0.46, 0.33))) return 0;
  223       if (!compare_point(f, f1, vec(1.0, 1.0))) return 0;
  224       const double new_energy = f.field_energy();
  225       if (!compare(new_energy, f1.field_energy(), "   total energy")) return 0;
  226       if (new_energy > last_energy * 1e-6) {
  227         master_printf("Energy decaying too slowly: from %g to %g (%g)\n", last_energy, new_energy,
  228                       new_energy / last_energy);
  229         return 0;
  230       } else {
  231         master_printf("Got newE/oldE of %g\n", new_energy / last_energy);
  232       }
  233       field_energy_check_time += deltaT;
  234     }
  235   }
  236   return 1;
  237 }
  238 
  239 int test_pml_tm(double eps(const vec &), int splitting, const char *mydirname) {
  240   double a = 10.0;
  241 
  242   grid_volume gv = voltwo(3.0, 3.0, a);
  243   structure s1(gv, eps, pml(1.0));
  244   structure s(gv, eps, pml(1.0), identity(), splitting);
  245   s.set_output_directory(mydirname);
  246   s1.set_output_directory(mydirname);
  247 
  248   master_printf("Testing TM pml while splitting into %d chunks...\n", splitting);
  249   fields f(&s);
  250   f.add_point_source(Ez, 0.8, 1.6, 0.0, 4.0, vec(1.299, 1.401), 1.0);
  251   fields f1(&s1);
  252   f1.add_point_source(Ez, 0.8, 1.6, 0.0, 4.0, vec(1.299, 1.401), 1.0);
  253   const double deltaT = 100.0;
  254   const double ttot = 3.1 * deltaT;
  255   double field_energy_check_time = deltaT;
  256 
  257   while (f.time() < f.last_source_time())
  258     f.step();
  259   while (f1.time() < f1.last_source_time())
  260     f1.step();
  261 
  262   double last_energy = f.field_energy();
  263   while (f.time() < ttot) {
  264     f.step();
  265     f1.step();
  266     if (f.time() >= field_energy_check_time) {
  267       if (!compare_point(f, f1, vec(0.5, 0.01))) return 0;
  268       if (!compare_point(f, f1, vec(0.46, 0.33))) return 0;
  269       if (!compare_point(f, f1, vec(1.0, 1.0))) return 0;
  270       const double new_energy = f.field_energy();
  271       if (!compare(new_energy, f1.field_energy(), "   total energy")) return 0;
  272       if (new_energy > last_energy * 4e-6) {
  273         master_printf("Energy decaying too slowly: from %g to %g (%g)\n", last_energy, new_energy,
  274                       new_energy / last_energy);
  275         return 0;
  276       } else {
  277         master_printf("Got newE/oldE of %g\n", new_energy / last_energy);
  278       }
  279       field_energy_check_time += deltaT;
  280     }
  281   }
  282   return 1;
  283 }
  284 
  285 int test_pml_te(double eps(const vec &), int splitting, const char *mydirname) {
  286   double a = 10.0;
  287 
  288   grid_volume gv = voltwo(3.0, 3.0, a);
  289   structure s1(gv, eps, pml(1.0));
  290   structure s(gv, eps, pml(1.0), identity(), splitting);
  291   s.set_output_directory(mydirname);
  292   s1.set_output_directory(mydirname);
  293 
  294   master_printf("Testing TE pml while splitting into %d chunks...\n", splitting);
  295   fields f(&s);
  296   f.add_point_source(Hz, 0.7, 1.5, 0.0, 4.0, vec(1.5, 1.5), 1.0);
  297   f.add_point_source(Hz, 0.7, 1.5, 0.0, 4.0, vec(1.37, 1.27), 1.0);
  298   fields f1(&s1);
  299   f1.add_point_source(Hz, 0.7, 1.5, 0.0, 4.0, vec(1.5, 1.5), 1.0);
  300   f1.add_point_source(Hz, 0.7, 1.5, 0.0, 4.0, vec(1.37, 1.27), 1.0);
  301   const double deltaT = 100.0;
  302   const double ttot = 3.1 * deltaT;
  303   double field_energy_check_time = deltaT;
  304 
  305   while (f.time() < f.last_source_time())
  306     f.step();
  307   while (f1.time() < f1.last_source_time())
  308     f1.step();
  309 
  310   double last_energy = f.field_energy();
  311   while (f.time() < ttot) {
  312     f.step();
  313     f1.step();
  314     if (f.time() >= field_energy_check_time) {
  315       if (!compare_point(f, f1, vec(0.5, 0.01))) return 0;
  316       if (!compare_point(f, f1, vec(0.46, 0.33))) return 0;
  317       if (!compare_point(f, f1, vec(1.0, 1.0))) return 0;
  318       const double new_energy = f.field_energy();
  319       if (!compare(new_energy, f1.field_energy(), "   total energy")) return 0;
  320       if (new_energy > last_energy * 1.1e-6) {
  321         master_printf("Energy decaying too slowly: from %g to %g (%g)\n", last_energy, new_energy,
  322                       new_energy / last_energy);
  323         return 0;
  324       } else {
  325         master_printf("Got newE/oldE of %g\n", new_energy / last_energy);
  326       }
  327       field_energy_check_time += deltaT;
  328     }
  329   }
  330   return 1;
  331 }
  332 
  333 int main(int argc, char **argv) {
  334   initialize mpi(argc, argv);
  335   quiet = true;
  336   const char *mydirname = "two_dimensional-out";
  337   trash_output_directory(mydirname);
  338   master_printf("Testing 2D...\n");
  339 
  340   for (int s = 2; s < 4; s++)
  341     if (!test_pml(one, s, mydirname)) abort("error in test_pml vacuum\n");
  342 
  343   for (int s = 2; s < 4; s++)
  344     if (!test_pml_tm(one, s, mydirname)) abort("error in test_pml_tm vacuum\n");
  345 
  346   for (int s = 2; s < 4; s++)
  347     if (!test_pml_te(one, s, mydirname)) abort("error in test_pml_te vacuum\n");
  348 
  349   for (int s = 2; s < 4; s++)
  350     if (!test_metal(one, s, mydirname)) abort("error in test_metal vacuum\n");
  351   // if (!test_metal(one, 200, mydirname)) abort("error in test_metal vacuum\n");
  352 
  353   for (int s = 2; s < 5; s++)
  354     if (!test_metal(targets, s, mydirname)) abort("error in test_metal targets\n");
  355   // if (!test_metal(targets, 60, mydirname)) abort("error in test_metal targets\n");
  356 
  357   for (int s = 2; s < 5; s++)
  358     if (!test_periodic(targets, s, mydirname)) abort("error in test_periodic targets\n");
  359   // if (!test_periodic(one, 200, mydirname))
  360   //  abort("error in test_periodic targets\n");
  361 
  362   for (int s = 2; s < 4; s++)
  363     if (!test_periodic_tm(one, s, mydirname)) abort("error in test_periodic_tm vacuum\n");
  364 
  365   return 0;
  366 }