"Fossies" - the Fresh Open Source Software Archive

Member "meep-1.10.0/src/structure_dump.cpp" (16 Mar 2019, 21952 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. For more information about "structure_dump.cpp" see the Fossies "Dox" file reference documentation and the last Fossies "Diffs" side-by-side code changes report: 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 // Dump/load raw structure data to/from an HDF5 file.  Only
   19 // works if the number of processors/chunks is the same.
   20 
   21 #include <stdio.h>
   22 #include <stdlib.h>
   23 #include <math.h>
   24 #include <string.h>
   25 
   26 #include "meep.hpp"
   27 #include "meep_internals.hpp"
   28 
   29 using namespace std;
   30 
   31 namespace meep {
   32 
   33 // Write the parameters required to reconstruct the susceptibility (id, noise_amp (for noisy),
   34 // omega_0, gamma, no_omega_0_denominator)
   35 void structure::write_susceptibility_params(h5file *file, const char *dname, int EorH) {
   36   // Get number of susceptibility params from first chunk, since all chunks will have
   37   // the same susceptibility list.
   38   size_t params_ntotal = 0;
   39   susceptibility *sus = chunks[0]->chiP[EorH];
   40   while (sus) {
   41     params_ntotal += sus->get_num_params() + 1;
   42     sus = sus->next;
   43   }
   44 
   45   // Write params
   46   size_t params_start = 0;
   47   file->create_data(dname, 1, &params_ntotal);
   48   if (am_master()) {
   49     susceptibility *sus = chunks[0]->chiP[EorH];
   50     while (sus) {
   51       sus->dump_params(file, &params_start);
   52       sus = sus->next;
   53     }
   54   }
   55 }
   56 
   57 void structure::dump_chunk_layout(const char *filename) {
   58   // Write grid_volume info for each chunk so we can reconstruct chunk division from split_by_cost
   59   size_t sz = num_chunks * 3;
   60   realnum *origins = new realnum[sz];
   61   size_t *nums = new size_t[sz];
   62   memset(nums, 0, sizeof(size_t) * sz);
   63   memset(origins, 0, sizeof(realnum) * sz);
   64   for (int i = 0; i < num_chunks; ++i) {
   65     int idx = i * 3;
   66     LOOP_OVER_DIRECTIONS(gv.dim, d) {
   67       origins[idx++] = chunks[i]->gv.origin_in_direction(d);
   68       nums[i * 3 + ((int)d % 3)] = chunks[i]->gv.num_direction(d);
   69     }
   70   }
   71   h5file file(filename, h5file::WRITE, true);
   72   file.create_data("gv_origins", 1, &sz);
   73   if (am_master()) {
   74     size_t gv_origins_start = 0;
   75     file.write_chunk(1, &gv_origins_start, &sz, origins);
   76   }
   77   file.create_data("gv_nums", 1, &sz);
   78   if (am_master()) {
   79     size_t nums_start = 0;
   80     file.write_chunk(1, &nums_start, &sz, nums);
   81   }
   82   delete[] origins;
   83   delete[] nums;
   84 }
   85 
   86 void structure::dump(const char *filename) {
   87   if (!quiet) master_printf("creating epsilon output file \"%s\"...\n", filename);
   88 
   89   // make/save a num_chunks x NUM_FIELD_COMPONENTS x 5 array counting
   90   // the number of entries in the chi1inv array for each chunk.
   91   size_t *num_chi1inv_ = new size_t[num_chunks * NUM_FIELD_COMPONENTS * 5];
   92   memset(num_chi1inv_, 0, sizeof(size_t) * size_t(num_chunks * NUM_FIELD_COMPONENTS * 5));
   93   size_t my_ntot = 0;
   94   for (int i = 0; i < num_chunks; i++)
   95     if (chunks[i]->is_mine()) {
   96       size_t ntot = chunks[i]->gv.ntot();
   97       for (int c = 0; c < NUM_FIELD_COMPONENTS; ++c)
   98         for (int d = 0; d < 5; ++d)
   99           if (chunks[i]->chi1inv[c][d])
  100             my_ntot += (num_chi1inv_[(i * NUM_FIELD_COMPONENTS + c) * 5 + d] = ntot);
  101     }
  102   size_t *num_chi1inv = new size_t[num_chunks * NUM_FIELD_COMPONENTS * 5];
  103   sum_to_master(num_chi1inv_, num_chi1inv, num_chunks * NUM_FIELD_COMPONENTS * 5);
  104   delete[] num_chi1inv_;
  105 
  106   // determine total dataset size and offset of this process's data
  107   size_t my_start = partial_sum_to_all(my_ntot) - my_ntot;
  108   size_t ntotal = sum_to_all(my_ntot);
  109 
  110   h5file file(filename, h5file::WRITE, true);
  111   size_t dims[3] = {(size_t)num_chunks, NUM_FIELD_COMPONENTS, 5};
  112   size_t start[3] = {0, 0, 0};
  113   file.create_data("num_chi1inv", 3, dims);
  114   if (am_master()) file.write_chunk(3, start, dims, num_chi1inv);
  115   delete[] num_chi1inv;
  116 
  117   // write the data
  118   file.create_data("chi1inv", 1, &ntotal);
  119   for (int i = 0; i < num_chunks; i++)
  120     if (chunks[i]->is_mine()) {
  121       size_t ntot = chunks[i]->gv.ntot();
  122       for (int c = 0; c < NUM_FIELD_COMPONENTS; ++c)
  123         for (int d = 0; d < 5; ++d)
  124           if (chunks[i]->chi1inv[c][d]) {
  125             file.write_chunk(1, &my_start, &ntot, chunks[i]->chi1inv[c][d]);
  126             my_start += ntot;
  127           }
  128     }
  129 
  130   // Get the sizes of susceptibility lists for chiP[E_stuff] and chiP[H_stuff]
  131   // Since this info is copied to each chunk, we can just get it from the first chunk
  132   size_t num_sus[2] = {0, 0};
  133   susceptibility *Esus = chunks[0]->chiP[E_stuff];
  134   while (Esus) {
  135     num_sus[E_stuff] += 1;
  136     Esus = Esus->next;
  137   }
  138   susceptibility *Hsus = chunks[0]->chiP[H_stuff];
  139   while (Hsus) {
  140     num_sus[H_stuff] += 1;
  141     Hsus = Hsus->next;
  142   }
  143 
  144   {
  145     // Write the number of susceptibilites
  146     size_t len = 2;
  147     file.create_data("num_sus", 1, &len);
  148     if (am_master()) {
  149       size_t start = 0;
  150       size_t ntot = 2;
  151       file.write_chunk(1, &start, &ntot, num_sus);
  152     }
  153   }
  154 
  155   // Get number of non-null sigma entries for each chiP in each chunk.
  156   // Assumes each susceptibility in the chiP[E_stuff] list has the
  157   // same number of non-null sigma elements. Likewise for chiP[H_stuff]
  158   size_t *my_num_sigmas[2];
  159   size_t *num_sigmas[2];
  160   for (int ft = 0; ft < 2; ++ft) {
  161     my_num_sigmas[ft] = new size_t[num_chunks];
  162     num_sigmas[ft] = new size_t[num_chunks];
  163     memset(my_num_sigmas[ft], 0, sizeof(size_t) * num_chunks);
  164     memset(num_sigmas[ft], 0, sizeof(size_t) * num_chunks);
  165   }
  166 
  167   for (int i = 0; i < num_chunks; ++i) {
  168     if (chunks[i]->is_mine()) {
  169       for (int ft = 0; ft < 2; ++ft) {
  170         susceptibility *sus = chunks[i]->chiP[ft];
  171         if (sus) {
  172           for (int c = 0; c < NUM_FIELD_COMPONENTS; ++c) {
  173             for (int d = 0; d < 5; ++d) {
  174               if (sus->sigma[c][d]) { my_num_sigmas[ft][i] += 1; }
  175             }
  176           }
  177         }
  178       }
  179     }
  180   }
  181 
  182   // Write num_sigmas data.
  183   {
  184     size_t ntot = num_chunks * 2;
  185     file.create_data("num_sigmas", 1, &ntot);
  186 
  187     for (int i = 0; i < num_chunks; ++i) {
  188       if (chunks[i]->is_mine()) {
  189         for (int ft = 0; ft < 2; ++ft) {
  190           size_t start = ft * num_chunks + i;
  191           size_t count = 1;
  192           file.write_chunk(1, &start, &count, &my_num_sigmas[ft][i]);
  193         }
  194       }
  195     }
  196   }
  197 
  198   file.prevent_deadlock();
  199   for (int ft = 0; ft < 2; ++ft) {
  200     sum_to_all(my_num_sigmas[ft], num_sigmas[ft], num_chunks);
  201   }
  202 
  203   size_t num_E_sigmas = 0;
  204   size_t num_H_sigmas = 0;
  205   for (int i = 0; i < num_chunks; ++i) {
  206     if (num_sigmas[E_stuff][i] != 0) { num_E_sigmas = num_sigmas[E_stuff][i]; }
  207     if (num_sigmas[H_stuff][i] != 0) { num_H_sigmas = num_sigmas[H_stuff][i]; }
  208   }
  209 
  210   // Allocate space for component and direction of non-null sigmas
  211   size_t *my_sigma_cd[2] = {NULL, NULL};
  212   my_sigma_cd[E_stuff] = new size_t[num_E_sigmas * 2];
  213   memset(my_sigma_cd[E_stuff], 0, sizeof(size_t) * num_E_sigmas * 2);
  214   my_sigma_cd[H_stuff] = new size_t[num_H_sigmas * 2];
  215   memset(my_sigma_cd[H_stuff], 0, sizeof(size_t) * num_H_sigmas * 2);
  216 
  217   size_t *sigma_cd[2] = {NULL, NULL};
  218   sigma_cd[E_stuff] = new size_t[num_E_sigmas * 2];
  219   memset(sigma_cd[E_stuff], 0, sizeof(size_t) * num_E_sigmas * 2);
  220   sigma_cd[H_stuff] = new size_t[num_H_sigmas * 2];
  221   memset(sigma_cd[H_stuff], 0, sizeof(size_t) * num_H_sigmas * 2);
  222 
  223   // Find component and direction of non-null sigmas
  224   {
  225     for (int ft = 0; ft < 2; ++ft) {
  226       int j = 0;
  227       bool done = false;
  228       for (int i = 0; i < num_chunks; ++i) {
  229         if (done) { break; }
  230         if (chunks[i]->is_mine()) {
  231           susceptibility *sus = chunks[i]->chiP[ft];
  232           if (sus) {
  233             for (int c = 0; c < NUM_FIELD_COMPONENTS; ++c) {
  234               for (int d = 0; d < 5; ++d) {
  235                 if (sus->sigma[c][d]) {
  236                   my_sigma_cd[ft][j] = c;
  237                   my_sigma_cd[ft][j + 1] = d;
  238                   j += 2;
  239                   done = true;
  240                 }
  241               }
  242             }
  243           }
  244         }
  245       }
  246     }
  247     bw_or_to_all(my_sigma_cd[E_stuff], sigma_cd[E_stuff], num_E_sigmas * 2);
  248     bw_or_to_all(my_sigma_cd[H_stuff], sigma_cd[H_stuff], num_H_sigmas * 2);
  249   }
  250 
  251   // Write location (component and direction) data of non-null sigmas (sigma[c][d])
  252   {
  253     size_t len = (num_E_sigmas + num_H_sigmas) * 2;
  254     file.create_data("sigma_cd", 1, &len);
  255     size_t start = 0;
  256     for (int ft = 0; ft < 2; ++ft) {
  257       if (am_master()) {
  258         size_t count = (ft == 0 ? num_E_sigmas : num_H_sigmas) * 2;
  259         file.write_chunk(1, &start, &count, sigma_cd[ft]);
  260         start += count;
  261       }
  262     }
  263   }
  264 
  265   // Write the actual data in a particular non-null sigma[c][d] for each susceptibility in this
  266   // chunk's chiP lists.
  267   size_t nsig[2] = {num_E_sigmas, num_H_sigmas};
  268   for (int i = 0; i < num_chunks; ++i) {
  269     for (int ft = 0; ft < 2; ++ft) {
  270       if (nsig[ft] != 0 && num_sigmas[ft][i]) {
  271         for (size_t j = 0; j < nsig[ft] * 2; j += 2) {
  272           char dname[20];
  273           int c = sigma_cd[ft][j];
  274           int d = sigma_cd[ft][j + 1];
  275           snprintf(dname, 20, "%c_%d_sigma_%d_%d", ft == 0 ? 'E' : 'H', i, c, d);
  276           size_t ntot = chunks[i]->gv.ntot() * num_sus[ft];
  277           file.create_data(dname, 1, &ntot);
  278           if (chunks[i]->is_mine()) {
  279             susceptibility *sus = chunks[i]->chiP[ft];
  280             size_t start = 0;
  281             while (sus) {
  282               size_t count = chunks[i]->gv.ntot();
  283               file.write_chunk(1, &start, &count, sus->sigma[c][d]);
  284               sus = sus->next;
  285               start += count;
  286             }
  287           }
  288         }
  289       }
  290     }
  291   }
  292 
  293   write_susceptibility_params(&file, "E_params", E_stuff);
  294   write_susceptibility_params(&file, "H_params", H_stuff);
  295 
  296   for (int ft = 0; ft < 2; ++ft) {
  297     delete[] sigma_cd[ft];
  298     delete[] my_sigma_cd[ft];
  299     delete[] num_sigmas[ft];
  300     delete[] my_num_sigmas[ft];
  301   }
  302 }
  303 
  304 // Reconstruct the chiP lists of susceptibilities from the params hdf5 data
  305 susceptibility *make_sus_list_from_params(h5file *file, int rank, size_t dims[3], size_t ntot) {
  306   susceptibility *sus = NULL;
  307   susceptibility *res = NULL;
  308   size_t start = 0;
  309 
  310   while (start < dims[0]) {
  311     size_t num_params_dims[3] = {1, 0, 0};
  312     realnum num_params;
  313     file->read_chunk(rank, &start, num_params_dims, &num_params);
  314     start += num_params_dims[0];
  315 
  316     if (num_params == 4) {
  317       // This is a lorentzian_susceptibility and the next 4 values in the dataset
  318       // are id, omega_0, gamma, and no_omega_0_denominator.
  319       size_t lorentzian_dims[3] = {4, 0, 0};
  320       realnum lorentzian_params[4];
  321       file->read_chunk(rank, &start, lorentzian_dims, lorentzian_params);
  322       start += lorentzian_dims[0];
  323 
  324       int id = (int)lorentzian_params[0];
  325       double omega_0 = lorentzian_params[1];
  326       double gamma = lorentzian_params[2];
  327       bool no_omega_0_denominator = (bool)lorentzian_params[3];
  328       if (sus) {
  329         sus->next = new lorentzian_susceptibility(omega_0, gamma, no_omega_0_denominator);
  330         sus->next->ntot = ntot;
  331         sus->next->set_id(id);
  332       } else {
  333         sus = new lorentzian_susceptibility(omega_0, gamma, no_omega_0_denominator);
  334         sus->ntot = ntot;
  335         sus->set_id(id);
  336         res = sus;
  337       }
  338       if (sus->next) sus = sus->next;
  339     } else if (num_params == 5) {
  340       // This is a noisy_lorentzian_susceptibility and the next 5 values in the dataset
  341       // are id, noise_amp, omega_0, gamma, and no_omega_0_denominator.
  342       size_t noisy_lorentzian_dims[3] = {5, 0, 0};
  343       realnum noisy_lorentzian_params[5];
  344       file->read_chunk(rank, &start, noisy_lorentzian_dims, noisy_lorentzian_params);
  345       start += noisy_lorentzian_dims[0];
  346 
  347       int id = (int)noisy_lorentzian_params[0];
  348       double noise_amp = noisy_lorentzian_params[1];
  349       double omega_0 = noisy_lorentzian_params[2];
  350       double gamma = noisy_lorentzian_params[3];
  351       bool no_omega_0_denominator = (bool)noisy_lorentzian_params[4];
  352       if (sus) {
  353         sus->next =
  354             new noisy_lorentzian_susceptibility(noise_amp, omega_0, gamma, no_omega_0_denominator);
  355         sus->next->ntot = ntot;
  356         sus->next->set_id(id);
  357       } else {
  358         sus =
  359             new noisy_lorentzian_susceptibility(noise_amp, omega_0, gamma, no_omega_0_denominator);
  360         sus->ntot = ntot;
  361         sus->set_id(id);
  362         res = sus;
  363       }
  364       if (sus->next) sus = sus->next;
  365     } else {
  366       abort("Invalid number of susceptibility parameters in structure::load");
  367     }
  368   }
  369   return res;
  370 }
  371 
  372 void structure::set_chiP_from_file(h5file *file, const char *dataset, field_type ft) {
  373   int rank = 0;
  374   size_t dims[3] = {0, 0, 0};
  375 
  376   file->read_size(dataset, &rank, dims, 1);
  377   if (rank != 1) abort("inconsistent data size in structure::load");
  378 
  379   if (dims[0] != 0) {
  380     for (int i = 0; i < num_chunks; ++i) {
  381       chunks[i]->chiP[ft] = make_sus_list_from_params(file, rank, dims, chunks[i]->gv.ntot());
  382     }
  383   }
  384 }
  385 
  386 void structure::load_chunk_layout(const char *filename, boundary_region &br) {
  387   // Load chunk grid_volumes from a file
  388   h5file file(filename, h5file::READONLY, true);
  389   size_t sz = num_chunks * 3;
  390   realnum *origins = new realnum[sz];
  391   memset(origins, 0, sizeof(realnum) * sz);
  392   size_t *nums = new size_t[sz];
  393   memset(nums, 0, sizeof(size_t) * sz);
  394 
  395   int origins_rank;
  396   size_t origins_dims;
  397   file.read_size("gv_origins", &origins_rank, &origins_dims, 1);
  398   if (origins_rank != 1 || origins_dims != sz) { abort("chunk mismatch in structure::load"); }
  399   if (am_master()) {
  400     size_t gv_origins_start = 0;
  401     file.read_chunk(1, &gv_origins_start, &origins_dims, origins);
  402   }
  403   file.prevent_deadlock();
  404   broadcast(0, origins, sz);
  405 
  406   int nums_rank;
  407   size_t nums_dims;
  408   file.read_size("gv_nums", &nums_rank, &nums_dims, 1);
  409   if (nums_rank != 1 || nums_dims != sz) { abort("chunk mismatch in structure::load"); }
  410   if (am_master()) {
  411     size_t gv_nums_start = 0;
  412     file.read_chunk(1, &gv_nums_start, &nums_dims, nums);
  413   }
  414   file.prevent_deadlock();
  415   broadcast(0, nums, sz);
  416 
  417   std::vector<grid_volume> gvs;
  418   // Populate a vector with the new grid_volumes
  419   for (int i = 0; i < num_chunks; ++i) {
  420     int idx = i * 3;
  421     grid_volume new_gv = gv;
  422     vec new_origin(new_gv.dim);
  423     LOOP_OVER_DIRECTIONS(gv.dim, d) {
  424       new_origin.set_direction(d, origins[idx++]);
  425       new_gv.set_num_direction(d, nums[i * 3 + ((int)d % 3)]);
  426     }
  427     new_gv.set_origin(new_origin);
  428     gvs.push_back(new_gv);
  429   }
  430 
  431   load_chunk_layout(gvs, br);
  432 
  433   delete[] origins;
  434   delete[] nums;
  435 }
  436 
  437 void structure::load_chunk_layout(const std::vector<grid_volume> &gvs, boundary_region &br) {
  438   // Recreate the chunks with the new grid_volumes
  439   for (int i = 0; i < num_chunks; ++i) {
  440     if (chunks[i]->refcount-- <= 1) delete chunks[i];
  441     int proc = i * count_processors() / num_chunks;
  442     chunks[i] = new structure_chunk(gvs[i], v, Courant, proc);
  443     br.apply(this, chunks[i]);
  444   }
  445   check_chunks();
  446 }
  447 
  448 void structure::load(const char *filename) {
  449   h5file file(filename, h5file::READONLY, true);
  450 
  451   if (!quiet) master_printf("reading epsilon from file \"%s\"...\n", filename);
  452 
  453   // make/save a num_chunks x NUM_FIELD_COMPONENTS x 5 array counting
  454   // the number of entries in the chi1inv array for each chunk.
  455   size_t *num_chi1inv = new size_t[num_chunks * NUM_FIELD_COMPONENTS * 5];
  456   int rank;
  457   size_t dims[3], _dims[3] = {(size_t)num_chunks, NUM_FIELD_COMPONENTS, 5};
  458   size_t start[3] = {0, 0, 0};
  459   file.read_size("num_chi1inv", &rank, dims, 3);
  460   if (rank != 3 || _dims[0] != dims[0] || _dims[1] != dims[1] || _dims[2] != dims[2])
  461     abort("chunk mismatch in structure::load");
  462   if (am_master()) file.read_chunk(3, start, dims, num_chi1inv);
  463 
  464   file.prevent_deadlock();
  465   broadcast(0, num_chi1inv, dims[0] * dims[1] * dims[2]);
  466 
  467   changing_chunks();
  468 
  469   // allocate data as needed and check sizes
  470   size_t my_ntot = 0;
  471   for (int i = 0; i < num_chunks; i++)
  472     if (chunks[i]->is_mine()) {
  473       size_t ntot = chunks[i]->gv.ntot();
  474       for (int c = 0; c < NUM_FIELD_COMPONENTS; ++c)
  475         for (int d = 0; d < 5; ++d) {
  476           size_t n = num_chi1inv[(i * NUM_FIELD_COMPONENTS + c) * 5 + d];
  477           if (n == 0) {
  478             delete[] chunks[i]->chi1inv[c][d];
  479             chunks[i]->chi1inv[c][d] = NULL;
  480           } else {
  481             if (n != ntot) abort("grid size mismatch %zd vs %zd in structure::load", n, ntot);
  482             chunks[i]->chi1inv[c][d] = new realnum[ntot];
  483             my_ntot += ntot;
  484           }
  485         }
  486     }
  487 
  488   // determine total dataset size and offset of this process's data
  489   size_t my_start = partial_sum_to_all(my_ntot) - my_ntot;
  490   size_t ntotal = sum_to_all(my_ntot);
  491 
  492   // read the data
  493   file.read_size("chi1inv", &rank, dims, 1);
  494   if (rank != 1 || dims[0] != ntotal) abort("inconsistent data size in structure::load");
  495   for (int i = 0; i < num_chunks; i++)
  496     if (chunks[i]->is_mine()) {
  497       size_t ntot = chunks[i]->gv.ntot();
  498       for (int c = 0; c < NUM_FIELD_COMPONENTS; ++c)
  499         for (int d = 0; d < 5; ++d)
  500           if (chunks[i]->chi1inv[c][d]) {
  501             file.read_chunk(1, &my_start, &ntot, chunks[i]->chi1inv[c][d]);
  502             my_start += ntot;
  503           }
  504     }
  505   // Create susceptibilites from params datasets
  506   set_chiP_from_file(&file, "E_params", E_stuff);
  507   set_chiP_from_file(&file, "H_params", H_stuff);
  508 
  509   // Read the number of susceptibilities in the chiP lists.
  510   size_t num_sus[] = {0, 0};
  511   {
  512     int rank = 0;
  513     size_t dims[] = {0, 0, 0};
  514     file.read_size("num_sus", &rank, dims, 1);
  515     if (dims[0] > 0) {
  516       if (am_master()) {
  517         size_t start = 0;
  518         size_t count = 2;
  519         file.read_chunk(rank, &start, &count, num_sus);
  520       }
  521     }
  522     file.prevent_deadlock();
  523     broadcast(0, num_sus, 2);
  524   }
  525 
  526   // Allocate and read non-null sigma entry data
  527   size_t *num_sigmas[2] = {NULL, NULL};
  528   num_sigmas[E_stuff] = new size_t[num_chunks];
  529   memset(num_sigmas[E_stuff], 0, sizeof(size_t) * num_chunks);
  530   num_sigmas[H_stuff] = new size_t[num_chunks];
  531   memset(num_sigmas[H_stuff], 0, sizeof(size_t) * num_chunks);
  532 
  533   // Read num_sigmas data
  534   {
  535     int rank = 0;
  536     size_t dims[] = {0, 0, 0};
  537     file.read_size("num_sigmas", &rank, dims, 1);
  538     if (dims[0] != (size_t)num_chunks * 2) { abort("inconsistent data size in structure::load"); }
  539     if (am_master()) {
  540       size_t start = 0;
  541       size_t count = num_chunks;
  542       file.read_chunk(rank, &start, &count, num_sigmas[E_stuff]);
  543       start += count;
  544       file.read_chunk(rank, &start, &count, num_sigmas[H_stuff]);
  545     }
  546     file.prevent_deadlock();
  547     broadcast(0, num_sigmas[E_stuff], num_chunks);
  548     broadcast(0, num_sigmas[H_stuff], num_chunks);
  549   }
  550 
  551   // Allocate space for component and direction data of the non-null susceptibilities
  552   size_t *sigma_cd[2] = {NULL, NULL};
  553 
  554   size_t nsig[2] = {0, 0};
  555   for (int ft = 0; ft < 2; ++ft) {
  556     for (int i = 0; i < num_chunks; ++i) {
  557       if (num_sigmas[ft][i] != 0) { nsig[ft] = num_sigmas[ft][i]; }
  558     }
  559   }
  560 
  561   sigma_cd[E_stuff] = new size_t[nsig[E_stuff] * 2];
  562   memset(sigma_cd[E_stuff], 0, sizeof(size_t) * nsig[E_stuff] * 2);
  563   sigma_cd[H_stuff] = new size_t[nsig[H_stuff] * 2];
  564   memset(sigma_cd[H_stuff], 0, sizeof(size_t) * nsig[H_stuff] * 2);
  565 
  566   // Read the component/direction data
  567   {
  568     int rank;
  569     size_t dims[] = {0, 0, 0};
  570     file.read_size("sigma_cd", &rank, dims, 1);
  571     if (dims[0] != 2 * (nsig[E_stuff] + nsig[H_stuff])) {
  572       abort("inconsistent data size in structure::load");
  573     }
  574 
  575     if (am_master()) {
  576       size_t start = 0;
  577       for (int ft = 0; ft < 2; ++ft) {
  578         size_t count = nsig[ft] * 2;
  579         file.read_chunk(rank, &start, &count, sigma_cd[ft]);
  580         start += count;
  581       }
  582     }
  583     file.prevent_deadlock();
  584     broadcast(0, sigma_cd[E_stuff], nsig[E_stuff] * 2);
  585     broadcast(0, sigma_cd[H_stuff], nsig[H_stuff] * 2);
  586   }
  587 
  588   for (int ft = 0; ft < 2; ++ft) {
  589     for (int i = 0; i < num_chunks; ++i) {
  590       for (int c = 0; c < NUM_FIELD_COMPONENTS; ++c) {
  591         for (int d = 0; d < 5; ++d) {
  592           susceptibility *sus = chunks[i]->chiP[ft];
  593           while (sus) {
  594             for (size_t j = 0; j < nsig[ft] * 2; j += 2) {
  595               int _c = sigma_cd[ft][j];
  596               int _d = sigma_cd[ft][j + 1];
  597               sus->trivial_sigma[_c][_d] = false;
  598             }
  599             sus = sus->next;
  600           }
  601         }
  602       }
  603     }
  604   }
  605 
  606   // Load sigma data
  607   for (int ft = 0; ft < 2; ++ft) {
  608     for (int i = 0; i < num_chunks; ++i) {
  609       for (int c = 0; c < NUM_FIELD_COMPONENTS; ++c) {
  610         for (int d = 0; d < 5; ++d) {
  611           char dname[20];
  612           snprintf(dname, 20, "%c_%d_sigma_%d_%d", ft == 0 ? 'E' : 'H', i, c, d);
  613           if (file.dataset_exists(dname)) {
  614             int rank;
  615             size_t dims[3] = {0, 0, 0};
  616             file.read_size(dname, &rank, dims, 1);
  617             if (num_sigmas[ft][i] && chunks[i]->is_mine()) {
  618               susceptibility *sus = chunks[i]->chiP[ft];
  619               size_t start = 0;
  620               while (sus) {
  621                 size_t count = chunks[i]->gv.ntot();
  622                 if (sus->sigma[c][d]) { delete[] sus->sigma[c][d]; }
  623                 sus->sigma[c][d] = new realnum[count];
  624                 sus->trivial_sigma[c][d] = false;
  625                 file.read_chunk(rank, &start, &count, sus->sigma[c][d]);
  626                 sus = sus->next;
  627                 start += count;
  628               }
  629             }
  630           }
  631         }
  632       }
  633     }
  634   }
  635 
  636   for (int ft = 0; ft < 2; ++ft) {
  637     delete[] num_sigmas[ft];
  638     delete[] sigma_cd[ft];
  639   }
  640 }
  641 } // namespace meep