Get source amplitude from hdf5 file (#388)

* Basic amplitude file func working

* Add linear_interpolate to meep.hpp to avoid redundancy

* libmeep doesn't need libmeepgeom

* Add namespace to linear_interpolate

* Use LDADD for all C++ tests

* Link libmeep.so with libctlgeom.so

* Remove superfluous changes

* Free memory after amp func is done with hdf5 data.

* Updates after review

* Create `add_volume_source` that takes array and size
* Error checking for real and imag dims sizes
* Print name of missing dataset on error
* Base linear interpolation on `where_` instead of geometry_lattice
* Don't need to link with libctlgeom

libmeepgeom/meepgeom.cpp

@@ -317,57 +317,6 @@ bool is_metal(meep::field_type ft, const material_type *material) {
   }
 }
 
-/* Linearly interpolate a given point in a 3d grid of data.  The point
-   coordinates should be in the range [0,1], or at the very least [-1,2]
-   ... anything outside [0,1] is *mirror* reflected into [0,1] */
-meep::realnum linear_interpolate(
-		     meep::realnum rx, meep::realnum ry, meep::realnum rz,
-		     meep::realnum *data, int nx, int ny, int nz, int stride)
-{
-     int x, y, z, x2, y2, z2;
-     meep::realnum dx, dy, dz;
-
-     /* mirror boundary conditions for r just beyond the boundary */
-     if (rx < 0.0) rx = -rx; else if (rx > 1.0) rx = 1.0 - rx;
-     if (ry < 0.0) ry = -ry; else if (ry > 1.0) ry = 1.0 - ry;
-     if (rz < 0.0) rz = -rz; else if (rz > 1.0) rz = 1.0 - rz;
-
-     /* get the point corresponding to r in the epsilon array grid: */
-     x = rx * nx; if (x == nx) --x;
-     y = ry * ny; if (y == ny) --y;
-     z = rz * nz; if (z == nz) --z;
-
-     /* get the difference between (x,y,z) and the actual point
-        ... we shift by 0.5 to center the data points in the pixels */
-     dx = rx * nx - x - 0.5;
-     dy = ry * ny - y - 0.5;
-     dz = rz * nz - z - 0.5;
-
-     /* get the other closest point in the grid, with mirror boundaries: */
-     x2 = (dx >= 0.0 ? x + 1 : x - 1);
-     if (x2 < 0) x2++; else if (x2 == nx) x2--;
-     y2 = (dy >= 0.0 ? y + 1 : y - 1);
-     if (y2 < 0) y2++; else if (y2 == ny) y2--;
-     z2 = (dz >= 0.0 ? z + 1 : z - 1);
-     if (z2 < 0) z2++; else if (z2 == nz) z2--;
-
-     /* take abs(d{xyz}) to get weights for {xyz} and {xyz}2: */
-     dx = fabs(dx);
-     dy = fabs(dy);
-     dz = fabs(dz);
-
-     /* define a macro to give us data(x,y,z) on the grid,
-	in row-major order (the order used by HDF5): */
-#define D(x,y,z) (data[(((x)*ny + (y))*nz + (z)) * stride])
-
-     return(((D(x,y,z)*(1.0-dx) + D(x2,y,z)*dx) * (1.0-dy) +
-	     (D(x,y2,z)*(1.0-dx) + D(x2,y2,z)*dx) * dy) * (1.0-dz) +
-	    ((D(x,y,z2)*(1.0-dx) + D(x2,y,z2)*dx) * (1.0-dy) +
-	     (D(x,y2,z2)*(1.0-dx) + D(x2,y2,z2)*dx) * dy) * dz);
-
-#undef D
-}
-
 // return material of the point p from the file (assumed already read)
 void epsilon_file_material(material_data *md, vector3 p)
 {
@@ -385,10 +334,10 @@ void epsilon_file_material(material_data *md, vector3 p)
   double rz = geometry_lattice.size.z == 0
     ? 0 : 0.5 + (p.z-geometry_center.z) / geometry_lattice.size.z;
   mm->epsilon_diag.x = mm->epsilon_diag.y = mm->epsilon_diag.z =
-    linear_interpolate(rx, ry, rz, md->epsilon_data,
-		       md->epsilon_dims[0],
-                       md->epsilon_dims[1],
-                       md->epsilon_dims[2], 1);
+    meep::linear_interpolate(rx, ry, rz, md->epsilon_data,
+                             md->epsilon_dims[0],
+                             md->epsilon_dims[1],
+                             md->epsilon_dims[2], 1);
   mm->epsilon_offdiag.x.re = mm->epsilon_offdiag.y.re = mm->epsilon_offdiag.z.re = 0;
 }
 

libmeepgeom/meepgeom.hpp

@@ -23,9 +23,6 @@
 #include <math.h>
 #include <vector>
 
-#include <meep.hpp>
-#include <ctlgeom.h>
-
 #include "meep.hpp"
 #include "material_data.hpp"
 
@@ -151,8 +148,6 @@ material_type make_file_material(const char *eps_input_file);
 vector3 vec_to_vector3(const meep::vec &pt);
 meep::vec vector3_to_vec(const vector3 v3);
 
-meep::realnum linear_interpolate(meep::realnum rx, meep::realnum ry, meep::realnum rz,
-                                 meep::realnum *data, int nx, int ny, int nz, int stride);
 void epsilon_file_material(material_data *md, vector3 p);
 bool susceptibility_equal(const susceptibility &s1, const susceptibility &s2);
 bool susceptibility_list_equal(const susceptibility_list &s1, const susceptibility_list &s2);

python/simulation.py

@@ -934,19 +934,28 @@ def add_source(self, src):
                     src.amp_func
                 )
         else:
-            if src.amp_func is None:
+            if src.amp_func_file:
                 self.fields.add_volume_source(
                     src.component,
                     src.src.swigobj,
                     where,
+                    src.amp_func_file,
+                    src.amp_func_dataset,
                     src.amplitude * 1.0,
                 )
-            else:
+            elif src.amp_func:
                 self.fields.add_volume_source(
                     src.component,
                     src.src.swigobj,
                     where,
                     src.amp_func,
+                    src.amplitude * 1.0,
+                )
+            else:
+                self.fields.add_volume_source(
+                    src.component,
+                    src.src.swigobj,
+                    where,
                     src.amplitude * 1.0
                 )
 

python/source.py

@@ -13,13 +13,16 @@ def check_positive(prop, val):
 
 class Source(object):
 
-    def __init__(self, src, component, center, size=Vector3(), amplitude=1.0, amp_func=None):
+    def __init__(self, src, component, center, size=Vector3(), amplitude=1.0, amp_func=None,
+                 amp_func_file='', amp_func_dataset=''):
         self.src = src
         self.component = component
         self.center = center
         self.size = size
         self.amplitude = complex(amplitude)
         self.amp_func = amp_func
+        self.amp_func_file = amp_func_file
+        self.amp_func_dataset = amp_func_dataset
 
 
 class SourceTime(object):

python/tests/source.py

@@ -1,13 +1,17 @@
 from __future__ import division
 
 import math
+import os
 import unittest
 
 import meep as mp
 from meep.geom import Cylinder, Vector3
 from meep.source import EigenModeSource, ContinuousSource, GaussianSource
 
 
+data_dir = os.path.join(os.path.realpath(os.path.dirname(__file__)), 'data')
+
+
 class TestEigenModeSource(unittest.TestCase):
 
     def test_eig_lattice_defaults(self):
@@ -135,5 +139,36 @@ def my_src_func(t):
 
         self.assertAlmostEqual(fp, -0.021997617628500023 + 0j)
 
+
+class TestAmpFileFunc(unittest.TestCase):
+
+    def init_and_run(self, file_func):
+        cell = mp.Vector3(1, 1)
+        resolution = 10
+        fcen = 0.8
+        df = 0.02
+
+        def ones(p):
+            return 1 + 1j
+
+        if file_func:
+            fname = os.path.join(data_dir, 'amp_func_file.h5')
+            dset = 'amp_data'
+            sources = [mp.Source(mp.ContinuousSource(fcen, fwidth=df), component=mp.Ez, center=mp.Vector3(),
+                                 amp_func_file=fname, amp_func_dataset=dset)]
+        else:
+            sources = [mp.Source(mp.ContinuousSource(fcen, fwidth=df), component=mp.Ez, center=mp.Vector3(),
+                                 amp_func=ones)]
+
+        sim = mp.Simulation(cell_size=cell, resolution=resolution, sources=sources)
+        sim.run(until=200)
+        return sim.get_field_point(mp.Ez, mp.Vector3())
+
+    def test_amp_file_func(self):
+        field_point_amp_file = self.init_and_run(file_func=True)
+        field_point_amp_func = self.init_and_run(file_func=False)
+        self.assertAlmostEqual(field_point_amp_file, field_point_amp_func)
+
+
 if __name__ == '__main__':
     unittest.main()

scheme/structure.cpp

@@ -193,57 +193,6 @@ static void read_epsilon_file(const char *eps_input_file)
   }
 }
 
-/* Linearly interpolate a given point in a 3d grid of data.  The point
-   coordinates should be in the range [0,1], or at the very least [-1,2]
-   ... anything outside [0,1] is *mirror* reflected into [0,1] */
-static meep::realnum linear_interpolate(
-		     meep::realnum rx, meep::realnum ry, meep::realnum rz,
-		     meep::realnum *data, int nx, int ny, int nz, int stride)
-{
-     int x, y, z, x2, y2, z2;
-     meep::realnum dx, dy, dz;
-
-     /* mirror boundary conditions for r just beyond the boundary */
-     if (rx < 0.0) rx = -rx; else if (rx > 1.0) rx = 1.0 - rx;
-     if (ry < 0.0) ry = -ry; else if (ry > 1.0) ry = 1.0 - ry;
-     if (rz < 0.0) rz = -rz; else if (rz > 1.0) rz = 1.0 - rz;
-
-     /* get the point corresponding to r in the epsilon array grid: */
-     x = rx * nx; if (x == nx) --x;
-     y = ry * ny; if (y == ny) --y;
-     z = rz * nz; if (z == nz) --z;
-
-     /* get the difference between (x,y,z) and the actual point
-        ... we shift by 0.5 to center the data points in the pixels */
-     dx = rx * nx - x - 0.5;
-     dy = ry * ny - y - 0.5;
-     dz = rz * nz - z - 0.5;
-
-     /* get the other closest point in the grid, with mirror boundaries: */
-     x2 = (dx >= 0.0 ? x + 1 : x - 1);
-     if (x2 < 0) x2++; else if (x2 == nx) x2--;
-     y2 = (dy >= 0.0 ? y + 1 : y - 1);
-     if (y2 < 0) y2++; else if (y2 == ny) y2--;
-     z2 = (dz >= 0.0 ? z + 1 : z - 1);
-     if (z2 < 0) z2++; else if (z2 == nz) z2--;
-
-     /* take abs(d{xyz}) to get weights for {xyz} and {xyz}2: */
-     dx = fabs(dx);
-     dy = fabs(dy);
-     dz = fabs(dz);
-
-     /* define a macro to give us data(x,y,z) on the grid,
-	in row-major order (the order used by HDF5): */
-#define D(x,y,z) (data[(((x)*ny + (y))*nz + (z)) * stride])
-
-     return(((D(x,y,z)*(1.0-dx) + D(x2,y,z)*dx) * (1.0-dy) +
-	     (D(x,y2,z)*(1.0-dx) + D(x2,y2,z)*dx) * dy) * (1.0-dz) +
-	    ((D(x,y,z2)*(1.0-dx) + D(x2,y,z2)*dx) * (1.0-dy) +
-	     (D(x,y2,z2)*(1.0-dx) + D(x2,y2,z2)*dx) * dy) * dz);
-
-#undef D
-}
-
 // return material of the point p from the file (assumed already read)
 static void epsilon_file_material(material_type &m, vector3 p)
 {
@@ -259,8 +208,8 @@ static void epsilon_file_material(material_type &m, vector3 p)
   double rz = geometry_lattice.size.z == 0
     ? 0 : 0.5 + (p.z-geometry_center.z) / geometry_lattice.size.z;
   mm->epsilon_diag.x = mm->epsilon_diag.y = mm->epsilon_diag.z =
-    linear_interpolate(rx, ry, rz, epsilon_data,
-		       epsilon_dims[0], epsilon_dims[1], epsilon_dims[2], 1);
+    meep::linear_interpolate(rx, ry, rz, epsilon_data,
+                             epsilon_dims[0], epsilon_dims[1], epsilon_dims[2], 1);
   mm->epsilon_offdiag.x = mm->epsilon_offdiag.y = mm->epsilon_offdiag.z = 0;
 }
 

src/fields.cpp

@@ -649,4 +649,54 @@ void fields::unset_solve_cw_omega() {
     chunks[i]->unset_solve_cw_omega();
 }
 
+/* Linearly interpolate a given point in a 3d grid of data.  The point
+   coordinates should be in the range [0,1], or at the very least [-1,2]
+   ... anything outside [0,1] is *mirror* reflected into [0,1] */
+realnum linear_interpolate(realnum rx, realnum ry, realnum rz, realnum *data,
+                           int nx, int ny, int nz, int stride) {
+
+  int x, y, z, x2, y2, z2;
+  realnum dx, dy, dz;
+
+  /* mirror boundary conditions for r just beyond the boundary */
+  if (rx < 0.0) rx = -rx; else if (rx > 1.0) rx = 1.0 - rx;
+  if (ry < 0.0) ry = -ry; else if (ry > 1.0) ry = 1.0 - ry;
+  if (rz < 0.0) rz = -rz; else if (rz > 1.0) rz = 1.0 - rz;
+
+  /* get the point corresponding to r in the epsilon array grid: */
+  x = rx * nx; if (x == nx) --x;
+  y = ry * ny; if (y == ny) --y;
+  z = rz * nz; if (z == nz) --z;
+
+  /* get the difference between (x,y,z) and the actual point
+     ... we shift by 0.5 to center the data points in the pixels */
+  dx = rx * nx - x - 0.5;
+  dy = ry * ny - y - 0.5;
+  dz = rz * nz - z - 0.5;
+
+  /* get the other closest point in the grid, with mirror boundaries: */
+  x2 = (dx >= 0.0 ? x + 1 : x - 1);
+  if (x2 < 0) x2++; else if (x2 == nx) x2--;
+  y2 = (dy >= 0.0 ? y + 1 : y - 1);
+  if (y2 < 0) y2++; else if (y2 == ny) y2--;
+  z2 = (dz >= 0.0 ? z + 1 : z - 1);
+  if (z2 < 0) z2++; else if (z2 == nz) z2--;
+
+  /* take abs(d{xyz}) to get weights for {xyz} and {xyz}2: */
+  dx = fabs(dx);
+  dy = fabs(dy);
+  dz = fabs(dz);
+
+  /* define a macro to give us data(x,y,z) on the grid,
+     in row-major order (the order used by HDF5): */
+#define D(x,y,z) (data[(((x)*ny + (y))*nz + (z)) * stride])
+
+  return(((D(x,y,z)*(1.0-dx) + D(x2,y,z)*dx) * (1.0-dy) +
+          (D(x,y2,z)*(1.0-dx) + D(x2,y2,z)*dx) * dy) * (1.0-dz) +
+         ((D(x,y,z2)*(1.0-dx) + D(x2,y,z2)*dx) * (1.0-dy) +
+          (D(x,y2,z2)*(1.0-dx) + D(x2,y2,z2)*dx) * dy) * dz);
+
+#undef D
+}
+
 } // namespace meep

src/h5file.cpp

@@ -332,8 +332,9 @@ realnum *h5file::read(const char *dataname,
 	close_data_id = false;
       }
       else {
-	CHECK(dataset_exists(file_id, dataname),
-	      "missing dataset in HDF5 file");
+        if (!dataset_exists(file_id, dataname)) {
+          abort("missing dataset in HDF5 file: %s", dataname);
+        }
 	data_id = H5Dopen(file_id, dataname);
       }
     }

src/meep.hpp

@@ -1425,6 +1425,12 @@ class fields {
                         int is_continuous = 0);
   void add_point_source(component c, const src_time &src,
                         const vec &, std::complex<double> amp = 1.0);
+  void add_volume_source(component c, const src_time &src, const volume &where_,
+                         std::complex<double> *arr, size_t dims[3],
+                         std::complex<double> amp);
+  void add_volume_source(component c, const src_time &src,
+                         const volume &where_, const char *filename,
+                         const char *dataset, std::complex<double> amp);
   void add_volume_source(component c, const src_time &src,
 			 const volume &,
 			 std::complex<double> A(const vec &),
@@ -1786,6 +1792,8 @@ std::complex<double> eigenmode_amplitude(void *vedata,
 double get_group_velocity(void *vedata);
 vec get_k(void *vedata);
 
+realnum linear_interpolate(realnum rx, realnum ry, realnum rz, realnum *data,
+                           int nx, int ny, int nz, int stride);
 } /* namespace meep */
 
 #endif /* MEEP_H */