"Fossies" - the Fresh Open Source Software Archive

Member "meep-1.10.0/tests/ring-ll.cpp" (11 Jan 2019, 6734 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.

    1 /***************************************************************/
    2 /* C++ port of meep/examples/ring.ctl, using the               */
    3 /* "low-level" meep C++ interface stack, which consists of     */
    4 /* libmeep_geom + libctlgeom + libmeep                         */
    5 /***************************************************************/
    6 /*
    7 ; Calculating 2d ring-resonator modes, from the Meep tutorial.
    8 */
    9 
   10 #include <stdio.h>
   11 #include <stdlib.h>
   12 #include <complex>
   13 #include <vector>
   14 
   15 #include "meep.hpp"
   16 
   17 #include "ctl-math.h"
   18 #include "ctlgeom.h"
   19 
   20 #include "meepgeom.hpp"
   21 
   22 using namespace meep;
   23 
   24 typedef std::complex<double> cdouble;
   25 
   26 vector3 v3(double x, double y = 0.0, double z = 0.0) {
   27   vector3 v;
   28   v.x = x;
   29   v.y = y;
   30   v.z = z;
   31   return v;
   32 }
   33 
   34 /***************************************************************/
   35 /* dummy material function needed to pass to structure( )      */
   36 /* constructor as a placeholder before we can call             */
   37 /* set_materials_from_geometry                                 */
   38 /***************************************************************/
   39 double dummy_eps(const vec &) { return 1.0; }
   40 
   41 /***************************************************************/
   42 /***************************************************************/
   43 /***************************************************************/
   44 int main(int argc, char *argv[]) {
   45   initialize mpi(argc, argv);
   46 
   47   double n = 3.4; // index of waveguide
   48   double w = 1.0; // width of waveguide
   49   double r = 1.0; // inner radius of ring
   50 
   51   double pad = 4;  // padding between waveguide and edge of PML
   52   double dpml = 2; // thickness of PML
   53 
   54   double sxy = 2.0 * (r + w + pad + dpml); // cell size
   55   double resolution = 10.0;
   56 
   57   // (set-param! resolution 10)
   58   // (set! geometry-lattice (make lattice (size sxy sxy no-size)))
   59   geometry_lattice.size.x = sxy;
   60   geometry_lattice.size.y = sxy;
   61   geometry_lattice.size.z = 0.0;
   62   grid_volume gv = voltwo(sxy, sxy, resolution);
   63   gv.center_origin();
   64 
   65   // (set! symmetries (list (make mirror-sym (direction Y))))
   66   // symmetry sym=mirror(Y, gv);
   67   symmetry sym = identity();
   68 
   69   // (set! pml-layers (list (make pml (thickness dpml))))
   70   // ; exploit the mirror symmetry in structure+source:
   71   structure the_structure(gv, dummy_eps, pml(dpml), sym);
   72 
   73   // ; Create a ring waveguide by two overlapping cylinders - later objects
   74   // ; take precedence over earlier objects, so we put the outer cylinder first.
   75   // ; and the inner (air) cylinder second.
   76   // (set! geometry (list
   77   //    (make cylinder (center 0 0) (height infinity)
   78   //        (radius (+ r w)) (material (make dielectric (index n))))
   79   //    (make cylinder (center 0 0) (height infinity)
   80   //        (radius r) (material air))))
   81   meep_geom::material_type dielectric = meep_geom::make_dielectric(n * n);
   82   geometric_object objects[2];
   83   vector3 v3zero = {0.0, 0.0, 0.0};
   84   vector3 zaxis = {0.0, 0.0, 1.0};
   85   objects[0] = make_cylinder(dielectric, v3zero, r + w, ENORMOUS, zaxis);
   86   objects[1] = make_cylinder(meep_geom::vacuum, v3zero, r, ENORMOUS, zaxis);
   87   geometric_object_list g = {2, objects};
   88   meep_geom::set_materials_from_geometry(&the_structure, g);
   89   fields f(&the_structure);
   90 
   91   // ; If we don't want to excite a specific mode symmetry, we can just
   92   // ; put a single point source at some arbitrary place, pointing in some
   93   // ; arbitrary direction.  We will only look for TM modes (E out of the plane).
   94   // (set! sources (list
   95   //              (make source
   96   //                (src (make gaussian-src (frequency fcen) (fwidth df)))
   97   //                (component Ez) (center (+ r 0.1) 0))))
   98   double fcen = 0.15; // ; pulse center frequency
   99   double df = 0.1;    // ; df
  100   gaussian_src_time src(fcen, df);
  101   volume v(vec(r + 0.1, 0.0), vec(0.0, 0.0));
  102   f.add_volume_source(Ez, src, v);
  103 
  104   // (run-sources+ 300
  105   //    (at-beginning output-epsilon)
  106   //    (after-sources (harminv Ez (vector3 (+ r 0.1)) fcen df)))
  107   while (f.round_time() < f.last_source_time())
  108     f.step();
  109 
  110   double T = 300.0;
  111   double stop_time = f.round_time() + T;
  112   std::vector<cdouble> fieldData;
  113   vec eval_pt(r + 0.1, 0.0);
  114   while (f.round_time() < stop_time) {
  115     f.step();
  116     fieldData.push_back(f.get_field(Ez, eval_pt));
  117   };
  118 
  119 #define MAXBANDS 100
  120   cdouble amp[MAXBANDS];
  121   double freq_re[MAXBANDS];
  122   double freq_im[MAXBANDS];
  123   double err[MAXBANDS];
  124   master_printf("starting do_harminv...\n");
  125   int bands = do_harminv(&fieldData[0], fieldData.size(), f.dt, fcen - 0.5 * df, fcen + 0.5 * df,
  126                          MAXBANDS, amp, freq_re, freq_im, err);
  127   master_printf("harminv0: | real(freq) | imag(freq)  |     Q      |  abs(amp)  |          amp     "
  128                 "     | err\n");
  129   master_printf("----------------------------------------------------------------------------------"
  130                 "----------\n");
  131   for (int nb = 0; nb < bands; nb++)
  132     master_printf("harminv0: | %.4e | %+.4e | %.4e | %.4e | {%+.2e,%+.2e} | %.1e \n", freq_re[nb],
  133                   freq_im[nb], 0.5 * freq_re[nb] / freq_im[nb], abs(amp[nb]), real(amp[nb]),
  134                   imag(amp[nb]), err[nb]);
  135 
  136   // test comparison with expected values
  137   int ref_bands = 3;
  138   double ref_freq_re[3] = {1.1807e-01, 1.4716e-01, 1.7525e-01};
  139   double ref_freq_im[3] = {-7.6133e-04, -2.1156e-04, -5.2215e-05};
  140   cdouble ref_amp[3] = {cdouble(-8.28e-04, -1.34e-03), cdouble(1.23e-03, -1.25e-02),
  141                         cdouble(2.83e-03, -6.52e-04)};
  142   if (bands != 3) abort("harminv found only %i/%i bands\n", bands, ref_bands);
  143   for (int nb = 0; nb < bands; nb++)
  144     if (fabs(freq_re[nb] - ref_freq_re[nb]) > 1.0e-2 * fabs(ref_freq_re[nb]) ||
  145         fabs(freq_im[nb] - ref_freq_im[nb]) > 1.0e-2 * fabs(ref_freq_im[nb]) ||
  146         abs(amp[nb] - ref_amp[nb]) > 1.0e-2 * abs(ref_amp[nb])
  147 
  148     )
  149       abort("harminv band %i disagrees with ref: {re f, im f, re A, im A}={%e,%e,%e,%e}!= "
  150             "{%e,%e,%e,%e}\n",
  151             nb, freq_re[nb], freq_im[nb], real(amp[nb]), imag(amp[nb]), ref_freq_re[nb],
  152             ref_freq_im[nb], real(ref_amp[nb]), imag(ref_amp[nb]));
  153 
  154   master_printf("all harminv results match reference values\n");
  155 
  156   // ; Output fields for one period at the end.  (If we output
  157   // ; at a single time, we might accidentally catch the Ez field
  158   // ; when it is almost zero and get a distorted view.)
  159   // (run-until (/ 1 fcen) (at-every (/ 1 fcen 20) output-efield-z))
  160   double DeltaT = 1.0 / (20 * fcen);
  161   double NextOutputTime = f.round_time() + DeltaT;
  162   while (f.round_time() < 1.0 / fcen) {
  163     f.step();
  164     if (f.round_time() >= NextOutputTime) {
  165       f.output_hdf5(Ez, f.total_volume());
  166       NextOutputTime += DeltaT;
  167     };
  168   };
  169 
  170   // this seems to be necessary to prevent failures
  171   all_wait();
  172 
  173   // success if we made it here
  174   return 0;
  175 }