"Fossies" - the Fresh Open Source Software Archive

Member "meep-1.10.0/tests/array-slice-ll.cpp" (5 Jun 2019, 8791 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 "array-slice-ll.cpp": 1.8.0_vs_1.9.0.

    1 /***************************************************************/
    2 /***************************************************************/
    3 /***************************************************************/
    4 #include <stdio.h>
    5 #include <stdlib.h>
    6 #include <string.h>
    7 #include <complex>
    8 
    9 #include "meep.hpp"
   10 
   11 #include "ctl-math.h"
   12 #include "ctlgeom.h"
   13 
   14 #include "meepgeom.hpp"
   15 
   16 #ifndef DATADIR
   17 #define DATADIR "./"
   18 #endif
   19 
   20 using namespace meep;
   21 
   22 typedef std::complex<double> cdouble;
   23 
   24 vector3 v3(double x, double y = 0.0, double z = 0.0) {
   25   vector3 v;
   26   v.x = x;
   27   v.y = y;
   28   v.z = z;
   29   return v;
   30 }
   31 
   32 // passthrough field function
   33 cdouble default_field_function(const cdouble *fields, const vec &loc, void *data_) {
   34   (void)loc;   // unused
   35   (void)data_; // unused
   36   return fields[0];
   37 }
   38 
   39 /***************************************************************/
   40 /***************************************************************/
   41 /***************************************************************/
   42 #define RELTOL 1.0e-6
   43 double Compare(double *d1, double *d2, int N, const char *Name) {
   44   double Norm1 = 0.0, Norm2 = 0.0, NormDelta = 0.0;
   45   for (int n = 0; n < N; n++) {
   46     Norm1 += d1[n] * d1[n];
   47     Norm2 += d2[n] * d2[n];
   48     NormDelta += (d1[n] - d2[n]) * (d1[n] - d2[n]);
   49   };
   50   Norm1 = sqrt(Norm1);
   51   Norm2 = sqrt(Norm2);
   52   NormDelta = sqrt(NormDelta);
   53   double RelErr = NormDelta / (0.5 * (Norm1 + Norm2));
   54   if (RelErr > RELTOL) abort("fail: rel error in %s data = %e\n", Name, RelErr);
   55   return RelErr;
   56 }
   57 
   58 double Compare(cdouble *d1, cdouble *d2, int N, const char *Name) {
   59   double Norm1 = 0.0, Norm2 = 0.0, NormDelta = 0.0;
   60   for (int n = 0; n < N; n++) {
   61     Norm1 += norm(d1[n]);
   62     Norm2 += norm(d2[n]);
   63     NormDelta += norm(d1[n] - d2[n]);
   64   };
   65   Norm1 = sqrt(Norm1);
   66   Norm2 = sqrt(Norm2);
   67   NormDelta = sqrt(NormDelta);
   68   double RelErr = NormDelta / (0.5 * (Norm1 + Norm2));
   69   if (RelErr > RELTOL) abort("fail: rel error in %s data = %e\n", Name, RelErr);
   70   return RelErr;
   71 }
   72 
   73 /***************************************************************/
   74 /* dummy material function needed to pass to structure( )      */
   75 /* constructor as a placeholder before we can call             */
   76 /* set_materials_from_geometry                                 */
   77 /***************************************************************/
   78 double dummy_eps(const vec &) { return 1.0; }
   79 
   80 /***************************************************************/
   81 /***************************************************************/
   82 /***************************************************************/
   83 void usage(char *progname) {
   84   master_printf("usage: %s [options]\n", progname);
   85   master_printf("options: \n");
   86   master_printf(" --use-symmetry    use geometric symmetries\n");
   87   master_printf(" --write-files     write reference data files\n");
   88   abort();
   89 }
   90 
   91 /***************************************************************/
   92 /***************************************************************/
   93 /***************************************************************/
   94 int main(int argc, char *argv[]) {
   95   initialize mpi(argc, argv);
   96 
   97   /***************************************************************/
   98   /* parse command-line options **********************************/
   99   /***************************************************************/
  100   bool use_symmetry = false;
  101   bool write_files = false;
  102   for (int narg = 1; narg < argc; narg++) {
  103     if (argv[narg] == 0) continue;
  104     if (!strcasecmp(argv[narg], "--use-symmetry")) {
  105       use_symmetry = true;
  106       master_printf("Using symmetry.\n");
  107     } else if (!strcasecmp(argv[narg], "--write-files")) {
  108       write_files = true;
  109       master_printf("writing HDF5 data files");
  110     } else {
  111       master_printf("unknown command-line option %s (aborting)", argv[narg]);
  112       usage(argv[0]);
  113     };
  114   };
  115 
  116   /***************************************************************/
  117   /* initialize geometry, similar to holey_wvg_cavity   **********/
  118   /***************************************************************/
  119   double eps = 13.0; // dielectric constant of waveguide
  120   double w = 1.2;    // width of waveguide
  121   double r = 0.36;   // radius of holes
  122   double d = 1.4;    // defect spacing (ordinary spacing = 1)
  123   int N = 3;         // number of holes on either side of defect
  124   double sy = 6.0;   // size of cell in y direction (perpendicular to wvg.)
  125   double pad = 2.0;  // padding between last hole and PML edge
  126   double dpml = 1.0; // PML thickness
  127   double sx = 2.0 * (pad + dpml + N) + d - 1.0; // size of cell in x dir
  128   double resolution = 20.0;
  129   geometry_lattice.size.x = sx;
  130   geometry_lattice.size.y = sy;
  131   geometry_lattice.size.z = 0.0;
  132   grid_volume gv = voltwo(sx, sy, resolution);
  133   gv.center_origin();
  134   symmetry sym = use_symmetry ? -mirror(Y, gv) : identity();
  135   structure the_structure(gv, dummy_eps, pml(dpml), sym);
  136   meep_geom::material_type vacuum = meep_geom::vacuum;
  137   meep_geom::material_type dielectric = meep_geom::make_dielectric(eps);
  138   geometric_object objects[7];
  139   vector3 origin = v3(0.0, 0.0, 0.0);
  140   vector3 xhat = v3(1.0, 0.0, 0.0);
  141   vector3 yhat = v3(0.0, 1.0, 0.0);
  142   vector3 zhat = v3(0.0, 0.0, 1.0);
  143   vector3 size = v3(ENORMOUS, w, ENORMOUS);
  144   double x0 = 0.5 * d;
  145   double deltax = 1.0;
  146   double height = ENORMOUS;
  147   objects[0] = make_block(dielectric, origin, xhat, yhat, zhat, size);
  148   int no = 1;
  149   for (int n = 0; n < N; n++) {
  150     vector3 center = v3(x0 + n * deltax, 0.0, 0.0);
  151     objects[no++] = make_cylinder(vacuum, center, r, height, zhat);
  152   };
  153   for (int n = 0; n < N; n++) {
  154     vector3 center = v3(-x0 - n * deltax, 0.0, 0.0);
  155     objects[no++] = make_cylinder(vacuum, center, r, height, zhat);
  156   };
  157   geometric_object_list g = {no, objects};
  158   meep_geom::set_materials_from_geometry(&the_structure, g);
  159   fields f(&the_structure);
  160 
  161   /***************************************************************/
  162   /* add source and timestep until source has finished (no later)*/
  163   /***************************************************************/
  164   double fcen = 0.25; // pulse center frequency
  165   double df = 0.2;    // pulse width (in frequency)
  166   gaussian_src_time src(fcen, df);
  167   component src_cmpt = Hz;
  168   f.add_point_source(src_cmpt, src, vec(0.0, 0.0));
  169   while (f.round_time() < f.last_source_time())
  170     f.step();
  171 
  172   /***************************************************************/
  173   /***************************************************************/
  174   /***************************************************************/
  175   double xMin = -0.25 * sx, xMax = +0.25 * sx;
  176   double yMin = -0.15 * sy, yMax = +0.15 * sy;
  177   volume v1d(vec(xMin, 0.0), vec(xMax, 0.0));
  178   volume v2d(vec(xMin, yMin), vec(xMax, yMax));
  179 
  180   int rank;
  181   size_t dims1D[1], dims2D[2];
  182   direction dirs1D[1], dirs2D[2];
  183   cdouble *file_slice1d = 0;
  184   double *file_slice2d = 0;
  185 
  186 #define H5FILENAME DATADIR "array-slice-ll-ref"
  187 #define NX 126
  188 #define NY 38
  189   if (write_files) {
  190     h5file *file = f.open_h5file(H5FILENAME);
  191     f.output_hdf5(Hz, v1d, file);
  192     f.output_hdf5(Sy, v2d, file);
  193     master_printf("Wrote binary data to file %s.h5\n", H5FILENAME);
  194     delete file;
  195     exit(0);
  196   } else {
  197     //
  198     // read 1D and 2D array-slice data from HDF5 file
  199     //
  200     h5file *file = f.open_h5file(H5FILENAME, h5file::READONLY);
  201     double *rdata = file->read("hz.r", &rank, dims1D, 1);
  202     if (rank != 1 || dims1D[0] != NX)
  203       abort("failed to read 1D data(hz.r) from file %s.h5", H5FILENAME);
  204     double *idata = file->read("hz.i", &rank, dims1D, 1);
  205     if (rank != 1 || dims1D[0] != NX)
  206       abort("failed to read 1D data(hz.i) from file %s.h5", H5FILENAME);
  207     file_slice1d = new cdouble[dims1D[0]];
  208     for (size_t n = 0; n < dims1D[0]; n++)
  209       file_slice1d[n] = cdouble(rdata[n], idata[n]);
  210     delete[] rdata;
  211     delete[] idata;
  212 
  213     file_slice2d = file->read("sy", &rank, dims2D, 2);
  214     if (rank != 2 || dims2D[0] != NX || dims2D[1] != NY)
  215       abort("failed to read 2D reference data from file %s.h5", H5FILENAME);
  216     delete file;
  217 
  218     //
  219     // generate 1D and 2D array slices and compare to
  220     // data read from file
  221     //
  222     rank = f.get_array_slice_dimensions(v1d, dims1D, dirs1D, true);
  223     if (rank != 1 || dims1D[0] != NX) abort("incorrect dimensions for 1D slice");
  224     cdouble *slice1d = f.get_complex_array_slice(v1d, Hz);
  225     double RelErr1D = Compare(slice1d, file_slice1d, NX, "Hz_1d");
  226     master_printf("1D: rel error %e\n", RelErr1D);
  227 
  228     rank = f.get_array_slice_dimensions(v2d, dims2D, dirs2D, true);
  229     if (rank != 2 || dims2D[0] != NX || dims2D[1] != NY) abort("incorrect dimensions for 2D slice");
  230     double *slice2d = f.get_array_slice(v2d, Sy);
  231     double RelErr2D = Compare(slice2d, file_slice2d, NX * NY, "Sy_2d");
  232     master_printf("2D: rel error %e\n", RelErr2D);
  233 
  234   }; // if (write_files) ... else ...
  235 
  236   return 0;
  237 }