"Fossies" - the Fresh Open Source Software Archive

Member "meep-1.10.0/tests/known_results.cpp" (27 Feb 2019, 6087 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 "known_results.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 
   21 #include <meep.hpp>
   22 using namespace meep;
   23 using std::complex;
   24 
   25 #include "config.h"
   26 
   27 double one(const vec &) { return 1.0; }
   28 
   29 double rods(const vec &r) {
   30   vec p = r;
   31   while (p.x() < -0.5)
   32     p.set_direction(X, p.x() + 1.0);
   33   while (p.x() > 0.5)
   34     p.set_direction(X, p.x() - 1.0);
   35   while (p.y() < -0.5)
   36     p.set_direction(Y, p.y() + 1.0);
   37   while (p.y() > 0.5)
   38     p.set_direction(Y, p.y() - 1.0);
   39   if (p.x() * p.x() + p.y() * p.y() < 0.3) return 12.0;
   40   return 1.0;
   41 }
   42 
   43 void compare(double b, double a, const char *n) {
   44   double thresh = sizeof(realnum) == sizeof(float) ? 1e-4 : 1e-5;
   45   if (fabs(a - b) > fabs(b) * thresh || b != b) {
   46     abort("Failed %s (%g instead of %g, relerr %0.2g)\n", n, a, b, fabs(a - b) / fabs(b));
   47   } else {
   48     master_printf("Passed %s\n", n);
   49   }
   50 }
   51 
   52 static double dpml = 1.0;
   53 
   54 double using_pml_ez(const grid_volume &gv, double eps(const vec &)) {
   55   const double ttot = 30.0;
   56   structure s(gv, eps, pml(dpml));
   57   fields f(&s);
   58   f.add_point_source(Ez, 0.2, 3.0, 0.0, 2.0, gv.center(), complex<double>(0, -2 * pi * 0.2));
   59   while (f.round_time() < ttot)
   60     f.step();
   61   monitor_point p;
   62   f.get_point(&p, gv.center());
   63   return real(p.get_component(Ez));
   64 }
   65 
   66 double x_periodic_y_pml(const grid_volume &gv, double eps(const vec &)) {
   67   const double ttot = 30.0;
   68   structure s(gv, eps, pml(dpml, Y));
   69   fields f(&s);
   70   f.add_point_source(Ez, 0.2, 3.0, 0.0, 2.0, gv.center(), complex<double>(0, -2 * pi * 0.2));
   71   f.use_bloch(X, 0.1);
   72   while (f.round_time() < ttot)
   73     f.step();
   74   monitor_point p;
   75   f.get_point(&p, gv.center());
   76   return real(p.get_component(Ez));
   77 }
   78 
   79 double x_periodic(const grid_volume &gv, double eps(const vec &)) {
   80   const double ttot = 30.0;
   81   structure s(gv, eps);
   82   fields f(&s);
   83   f.add_point_source(Ez, 0.2, 3.0, 0.0, 2.0, gv.center(), complex<double>(0, -2 * pi * 0.2));
   84   f.use_bloch(X, 0.1);
   85   while (f.round_time() < ttot)
   86     f.step();
   87   monitor_point p;
   88   f.get_point(&p, gv.center());
   89   return real(p.get_component(Ez));
   90 }
   91 
   92 double periodic_ez(const grid_volume &gv, double eps(const vec &)) {
   93   const double ttot = 30.0;
   94   structure s(gv, eps);
   95   fields f(&s);
   96   f.add_point_source(Ez, 0.2, 3.0, 0.0, 2.0, gv.center(), complex<double>(0, -2 * pi * 0.2));
   97   vec k;
   98   switch (gv.dim) {
   99     case D1: k = vec(0.3); break;
  100     case D2: k = vec(0.3, 0.4); break;
  101     case D3: k = vec(0.3, 0.5, 0.8); break;
  102     case Dcyl: k = veccyl(0.3, 0.2); break;
  103   }
  104   f.use_bloch(k);
  105   while (f.round_time() < ttot)
  106     f.step();
  107   monitor_point p;
  108   f.get_point(&p, gv.center());
  109   return real(p.get_component(Ez));
  110 }
  111 
  112 double metallic_ez(const grid_volume &gv, double eps(const vec &)) {
  113   const double ttot = 10.0;
  114   structure s(gv, eps);
  115   fields f(&s);
  116   f.add_point_source(Ez, 0.2, 3.0, 0.0, 2.0, gv.center(), complex<double>(0, -2 * pi * 0.2));
  117   while (f.round_time() < ttot)
  118     f.step();
  119   monitor_point p;
  120   f.get_point(&p, gv.center());
  121   return real(p.get_component(Ez));
  122 }
  123 
  124 double sigma(const vec &) { return 7.63; }
  125 
  126 double polariton_ex(const grid_volume &gv, double eps(const vec &)) {
  127   const double ttot = 10.0;
  128   structure s(gv, eps);
  129   s.add_susceptibility(sigma, E_stuff, lorentzian_susceptibility(0.3, 0.1));
  130   fields f(&s);
  131   f.add_point_source(Ex, 0.2, 3.0, 0.0, 2.0, gv.center(), complex<double>(0, -2 * pi * 0.2));
  132   while (f.round_time() < ttot)
  133     f.step();
  134   monitor_point p;
  135   f.get_point(&p, gv.center());
  136   return real(p.get_component(Ex));
  137 }
  138 
  139 double polariton_energy(const grid_volume &gv, double eps(const vec &)) {
  140   const double ttot = 10.0;
  141   structure s(gv, eps);
  142   s.add_susceptibility(sigma, E_stuff, lorentzian_susceptibility(0.3, 0.1));
  143   fields f(&s, 0);
  144   f.add_point_source(Ex, 0.2, 3.0, 0.0, 2.0, gv.center(), complex<double>(0, -2 * pi * 0.2));
  145   while (f.round_time() < ttot)
  146     f.step();
  147   return f.field_energy();
  148 }
  149 
  150 int main(int argc, char **argv) {
  151   initialize mpi(argc, argv);
  152   quiet = true;
  153   const char *mydirname = "known_results-out";
  154   trash_output_directory(mydirname);
  155   master_printf("Testing with some known results...\n");
  156   const double a = 10.0;
  157 
  158   compare(-0.0894851, polariton_ex(volone(1.0, a), one), "1D polariton");
  159   compare(0.0863443, polariton_energy(volone(1.0, a), one), "1D polariton energy");
  160   compare(5.20605, metallic_ez(voltwo(1.0, 1.0, a), one), "1x1 metallic 2D TM");
  161   compare(0.883776, using_pml_ez(voltwo(1.0 + 2 * dpml, 1.0 + 2 * dpml, a), one), "1x1 PML 2D TM");
  162   compare(0.110425, x_periodic(voltwo(1.0, 1.0, a), one), "1x1 X periodic 2D TM");
  163   compare(-4.78767, periodic_ez(voltwo(1.0, 3.0, a), rods), "1x1 fully periodic 2D TM rods");
  164   compare(1.12502, periodic_ez(voltwo(1.0, 3.0, a), one), "1x1 fully periodic 2D TM");
  165   compare(0.608815, x_periodic_y_pml(voltwo(1.0, 1.0 + 2 * dpml, a), one),
  166           "1x1 X periodic Y PML 2D TM");
  167   compare(-41.8057, metallic_ez(vol3d(1.0, 1.0, 1.0, a), one), "1x1x1 metallic 3D");
  168   compare(-100.758, x_periodic(vol3d(1.0, 1.0, 1.0, a), one), "1x1x1 X periodic 3D");
  169   compare(-101.398, x_periodic_y_pml(vol3d(1.0, 1.0 + 2 * dpml, 1.0, a), one),
  170           "1x1x1 X periodic Y PML 3D");
  171   compare(-103.844, periodic_ez(vol3d(1.0, 1.0, 1.0, a), rods), "1x1x1 fully periodic 3D rods");
  172   compare(-99.1618, periodic_ez(vol3d(1.0, 1.0, 1.0, a), one), "1x1x1 fully periodic 3D");
  173 
  174   return 0;
  175 }