From b5055af9d9228145f544df473467a05ccf7ef8b1 Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Thu, 17 May 2018 11:22:25 -0500 Subject: [PATCH 1/7] Add python Prism wrapper --- python/geom.py | 9 ++++++ python/meep.i | 14 +-------- python/tests/bend_flux.py | 22 +++++++++----- python/tests/geom.py | 9 ++++++ python/typemap_utils.cpp | 61 +++++++++++++++++++++++++++++++++++++++ 5 files changed, 95 insertions(+), 20 deletions(-) diff --git a/python/geom.py b/python/geom.py index 5a207bc5a..01ff7ff06 100644 --- a/python/geom.py +++ b/python/geom.py @@ -341,6 +341,15 @@ def __init__(self, **kwargs): super(Ellipsoid, self).__init__(**kwargs) +class Prism(GeometricObject): + + def __init__(self, vertices, height, axis=Vector3(z=1), **kwargs): + self.vertices = vertices + self.height = height + self.axis = axis + super(Prism, self).__init__(**kwargs) + + class Matrix(object): def __init__(self, c1=Vector3(), c2=Vector3(), c3=Vector3()): diff --git a/python/meep.i b/python/meep.i index f0898f0e7..929d56a9b 100644 --- a/python/meep.i +++ b/python/meep.i @@ -58,19 +58,6 @@ typedef struct { #include "typemap_utils.cpp" -static int get_attr_int(PyObject *py_obj, int *result, const char *name) { - PyObject *py_attr = PyObject_GetAttrString(py_obj, name); - - if (!py_attr) { - PyErr_Format(PyExc_ValueError, "Class attribute '%s' is None\n", name); - return 0; - } - - *result = PyInteger_AsLong(py_attr); - Py_DECREF(py_attr); - return 1; -} - static PyObject *py_source_time_object() { static PyObject *source_time_object = NULL; if (source_time_object == NULL) { @@ -1044,6 +1031,7 @@ void display_geometric_object_info(int indentby, GEOMETRIC_OBJECT o); Medium, NoisyDrudeSusceptibility, NoisyLorentzianSusceptibility, + Prism, Sphere, Susceptibility, Vector3, diff --git a/python/tests/bend_flux.py b/python/tests/bend_flux.py index e1c019409..5378891c5 100644 --- a/python/tests/bend_flux.py +++ b/python/tests/bend_flux.py @@ -17,13 +17,21 @@ def init(self, no_bend=False): wvg_xcen = 0.5 * (sx - w - (2 * pad)) if no_bend: - geometry = [mp.Block(mp.Vector3(mp.inf, w, mp.inf), center=mp.Vector3(0, wvg_ycen), - material=mp.Medium(epsilon=12))] + no_bend_vertices = [mp.Vector3(-0.5 * sx, wvg_ycen - 0.5 * w), + mp.Vector3(+0.5 * sx, wvg_ycen - 0.5 * w), + mp.Vector3(+0.5 * sx, wvg_ycen + 0.5 * w), + mp.Vector3(-0.5 * sx, wvg_ycen + 0.5 * w)] + + geometry = [mp.Prism(no_bend_vertices, 0, material=mp.Medium(epsilon=12))] else: - geometry = [mp.Block(mp.Vector3(sx - pad, w, mp.inf), center=mp.Vector3(-0.5 * pad, wvg_ycen), - material=mp.Medium(epsilon=12)), - mp.Block(mp.Vector3(w, sy - pad, mp.inf), center=mp.Vector3(wvg_xcen, 0.5 * pad), - material=mp.Medium(epsilon=12))] + bend_vertices = [mp.Vector3(-0.5 * sx, wvg_ycen - 0.5 * w), + mp.Vector3(wvg_xcen + 0.5 * w, wvg_ycen - 0.5 * w), + mp.Vector3(wvg_xcen + 0.5 * w, 0.5 * sy), + mp.Vector3(wvg_xcen - 0.5 * w, 0.5 * sy), + mp.Vector3(wvg_xcen - 0.5 * w, wvg_ycen + 0.5 * w), + mp.Vector3(-0.5 * sx, wvg_ycen + 0.5 * w)] + + geometry = [mp.Prism(bend_vertices, 0, material=mp.Medium(epsilon=12))] fcen = 0.15 df = 0.1 @@ -89,7 +97,7 @@ def test_bend_flux(self): res = list(zip(mp.get_flux_freqs(self.trans), mp.get_fluxes(self.trans), mp.get_fluxes(self.refl))) - np.testing.assert_allclose(expected, res[:20]) + np.testing.assert_allclose(expected, res[:20], rtol=1e-2) # Real run self.sim = None diff --git a/python/tests/geom.py b/python/tests/geom.py index e23623bcd..86d881839 100644 --- a/python/tests/geom.py +++ b/python/tests/geom.py @@ -245,6 +245,15 @@ def test_contains_point(self): self.assertIn(zeros(), e) +class TestPrism(unittest.TestCase): + + def test_contains_point(self): + vertices = [gm.Vector3(-1, 1), gm.Vector3(1, 1), gm.Vector3(1, -1), gm.Vector3(-1, -1)] + p = gm.Prism(vertices, height=1) + self.assertIn(zeros(), p) + self.assertNotIn(gm.Vector3(2, 2), p) + + class TestMedium(unittest.TestCase): def test_D_conductivity(self): diff --git a/python/typemap_utils.cpp b/python/typemap_utils.cpp index 646b098a0..6f219274b 100644 --- a/python/typemap_utils.cpp +++ b/python/typemap_utils.cpp @@ -263,6 +263,19 @@ static int get_attr_dbl(PyObject *py_obj, double *result, const char *name) { return 1; } +static int get_attr_int(PyObject *py_obj, int *result, const char *name) { + PyObject *py_attr = PyObject_GetAttrString(py_obj, name); + + if (!py_attr) { + PyErr_Format(PyExc_ValueError, "Class attribute '%s' is None\n", name); + return 0; + } + + *result = PyInteger_AsLong(py_attr); + Py_XDECREF(py_attr); + return 1; +} + static int get_attr_material(PyObject *po, material_type *m) { PyObject *py_material = PyObject_GetAttrString(po, "material"); @@ -541,6 +554,51 @@ static int pyellipsoid_to_ellipsoid(PyObject *py_ell, geometric_object *e) { return 1; } +static int pyprism_to_prism(PyObject *py_prism, geometric_object *p) { + material_type material; + int height; + vector3 axis; + + if (!get_attr_material(py_prism, &material) || + !get_attr_int(py_prism, &height, "height") || + !get_attr_v3(py_prism, &axis, "axis")) { + + return 0; + } + + PyObject *py_vert_list = PyObject_GetAttrString(py_prism, "vertices"); + + if (!py_vert_list) { + PyErr_PrintEx(0); + return 0; + } + + if (!PyList_Check(py_vert_list)) { + PyErr_SetString(PyExc_TypeError, "Expected Prism.vertices to be a list\n"); + return 0; + } + + int num_vertices = PyList_Size(py_vert_list); + vector3 *vertices = new vector3[num_vertices]; + + for (Py_ssize_t i = 0; i < num_vertices; ++i) { + vector3 v3; + if (!pyv3_to_v3(PyList_GetItem(py_vert_list, i), &v3)) { + return 0; + } + vertices[i].x = v3.x; + vertices[i].y = v3.y; + vertices[i].z = v3.z; + } + + *p = make_prism(material, vertices, num_vertices, height, axis); + + delete [] vertices; + Py_DECREF(py_vert_list); + + return 1; +} + static int py_gobj_to_gobj(PyObject *po, geometric_object *o) { int success = 0; std::string go_type = py_class_name_as_string(po); @@ -563,6 +621,9 @@ static int py_gobj_to_gobj(PyObject *po, geometric_object *o) { else if (go_type == "Ellipsoid") { success = pyellipsoid_to_ellipsoid(po, o); } + else if (go_type == "Prism") { + success = pyprism_to_prism(po, o); + } else { PyErr_Format(PyExc_TypeError, "Error: %s is not a valid GeometricObject type\n", go_type.c_str()); return 0; From bf8409f8d025683c0a637f93f307934768f87dac Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Thu, 17 May 2018 15:20:13 -0500 Subject: [PATCH 2/7] Increase tolerance and reuse compare_arrays from tests/mpb.py --- python/tests/bend_flux.py | 6 +++--- python/tests/mpb.py | 34 +++++++++++----------------------- python/tests/utils.py | 15 +++++++++++++++ 3 files changed, 29 insertions(+), 26 deletions(-) create mode 100644 python/tests/utils.py diff --git a/python/tests/bend_flux.py b/python/tests/bend_flux.py index 5378891c5..758ce3358 100644 --- a/python/tests/bend_flux.py +++ b/python/tests/bend_flux.py @@ -3,6 +3,7 @@ import unittest import numpy as np import meep as mp +from utils import compare_arrays class TestBendFlux(unittest.TestCase): @@ -32,7 +33,6 @@ def init(self, no_bend=False): mp.Vector3(-0.5 * sx, wvg_ycen + 0.5 * w)] geometry = [mp.Prism(bend_vertices, 0, material=mp.Medium(epsilon=12))] - fcen = 0.15 df = 0.1 sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), component=mp.Ez, @@ -97,7 +97,7 @@ def test_bend_flux(self): res = list(zip(mp.get_flux_freqs(self.trans), mp.get_fluxes(self.trans), mp.get_fluxes(self.refl))) - np.testing.assert_allclose(expected, res[:20], rtol=1e-2) + compare_arrays(self, np.array(expected), np.array(res[:20]), tol=1e-3) # Real run self.sim = None @@ -132,7 +132,7 @@ def test_bend_flux(self): res = list(zip(mp.get_flux_freqs(self.trans), mp.get_fluxes(self.trans), mp.get_fluxes(self.refl))) - np.testing.assert_allclose(expected, res[:20]) + compare_arrays(self, np.array(expected), np.array(res[:20]), tol=1e-2) if __name__ == '__main__': diff --git a/python/tests/mpb.py b/python/tests/mpb.py index 368a36096..34a5541d7 100644 --- a/python/tests/mpb.py +++ b/python/tests/mpb.py @@ -14,6 +14,7 @@ from scipy.optimize import ridder import meep as mp from meep import mpb +from utils import compare_arrays class TestModeSolver(unittest.TestCase): @@ -178,20 +179,7 @@ def check_fields_against_h5(self, ref_path, field, suffix=''): ref_arr[1::3] = ref_y.ravel() ref_arr[2::3] = ref_z.ravel() - self.compare_arrays(ref_arr, field) - - def compare_arrays(self, exp, res, tol=1e-3): - exp_1d = exp.ravel() - res_1d = res.ravel() - - norm_exp = np.linalg.norm(exp_1d) - norm_res = np.linalg.norm(res_1d) - - if norm_exp == 0: - self.assertEqual(norm_res, 0) - else: - diff = np.linalg.norm(res_1d - exp_1d) / norm_exp - self.assertLess(diff, tol) + compare_arrays(self, ref_arr, field) def compare_h5_files(self, ref_path, res_path, tol=1e-3): with h5py.File(ref_path) as ref: @@ -200,7 +188,7 @@ def compare_h5_files(self, ref_path, res_path, tol=1e-3): if k == 'description': self.assertEqual(ref[k].value, res[k].value) else: - self.compare_arrays(ref[k].value, res[k].value, tol=tol) + compare_arrays(self, ref[k].value, res[k].value, tol=tol) def test_update_band_range_data(self): brd = [] @@ -407,7 +395,7 @@ def test_compute_field_energy(self): expected_energy_in_dielectric = 0.6990769686037558 - self.compare_arrays(np.array(expected_energy), np.array(energy)) + compare_arrays(self, np.array(expected_energy), np.array(energy)) self.assertAlmostEqual(expected_energy_in_dielectric, energy_in_dielectric, places=3) def test_compute_group_velocity(self): @@ -693,7 +681,7 @@ def get_dpwr(ms, band): with h5py.File(ref_path, 'r') as f: expected = f['data-new'].value - self.compare_arrays(expected, converted_dpwr[-1]) + compare_arrays(self, expected, converted_dpwr[-1]) def test_hole_slab(self): from mpb_hole_slab import ms @@ -958,7 +946,7 @@ def test_tri_rods(self): expected_re = f['z.r-new'].value expected_im = f['z.i-new'].value expected = np.vectorize(complex)(expected_re, expected_im) - self.compare_arrays(expected, new_efield) + compare_arrays(self, expected, new_efield) ms.run_te() @@ -992,7 +980,7 @@ def test_tri_rods(self): with h5py.File(ref_path, 'r') as f: ref = f['data-new'].value - self.compare_arrays(ref, new_eps, tol=1e-3) + compare_arrays(self, ref, new_eps, tol=1e-3) def test_subpixel_averaging(self): ms = self.init_solver() @@ -1065,7 +1053,7 @@ def test_run_te_with_mu_material(self): mu = ms.get_mu() with h5py.File(data_path, 'r') as f: data = f['data'].value - self.compare_arrays(data, mu) + compare_arrays(self, data, mu) def test_output_tot_pwr(self): ms = self.init_solver() @@ -1083,7 +1071,7 @@ def test_output_tot_pwr(self): with h5py.File(ref_path, 'r') as f: expected = f['data'].value - self.compare_arrays(expected, arr) + compare_arrays(self, expected, arr) def test_get_eigenvectors(self): ms = self.init_solver() @@ -1196,7 +1184,7 @@ def test_epsilon_input_file(self): ] self.check_band_range_data(expected_brd, ms.band_range_data) - self.compare_arrays(expected_freqs, ms.all_freqs[-1]) + compare_arrays(self, expected_freqs, ms.all_freqs[-1]) self.check_gap_list(expected_gap_list, ms.gap_list) @@ -1340,7 +1328,7 @@ def test_multiply_bloch_in_map_data(self): md = mpb.MPBData(rectify=True, resolution=32, periods=3) result2 = md.convert(efield_no_bloch) - self.compare_arrays(result1, result2, tol=1e-7) + compare_arrays(self, result1, result2, tol=1e-7) def test_poynting(self): ms = self.init_solver() diff --git a/python/tests/utils.py b/python/tests/utils.py new file mode 100644 index 000000000..2e011d1f7 --- /dev/null +++ b/python/tests/utils.py @@ -0,0 +1,15 @@ +import numpy as np + + +def compare_arrays(test_instance, exp, res, tol=1e-3): + exp_1d = exp.ravel() + res_1d = res.ravel() + + norm_exp = np.linalg.norm(exp_1d) + norm_res = np.linalg.norm(res_1d) + + if norm_exp == 0: + test_instance.assertEqual(norm_res, 0) + else: + diff = np.linalg.norm(res_1d - exp_1d) / norm_exp + test_instance.assertLess(diff, tol) From f4b705666619f81910ebc990ec2df771542ad0f0 Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Fri, 18 May 2018 08:23:21 -0500 Subject: [PATCH 3/7] Simplify vector assignment --- python/typemap_utils.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/python/typemap_utils.cpp b/python/typemap_utils.cpp index 6f219274b..fd695fc43 100644 --- a/python/typemap_utils.cpp +++ b/python/typemap_utils.cpp @@ -586,9 +586,7 @@ static int pyprism_to_prism(PyObject *py_prism, geometric_object *p) { if (!pyv3_to_v3(PyList_GetItem(py_vert_list, i), &v3)) { return 0; } - vertices[i].x = v3.x; - vertices[i].y = v3.y; - vertices[i].z = v3.z; + vertices[i] = v3; } *p = make_prism(material, vertices, num_vertices, height, axis); From f425364afdeba1727a1a8aad4375439a9979c15d Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Wed, 23 May 2018 10:16:32 -0500 Subject: [PATCH 4/7] Get new scheme data with fixed size blocks --- python/tests/bend_flux.py | 88 +++++++++++++++++------------------ scheme/examples/bend-flux.ctl | 20 ++++---- 2 files changed, 53 insertions(+), 55 deletions(-) diff --git a/python/tests/bend_flux.py b/python/tests/bend_flux.py index 758ce3358..a5272073d 100644 --- a/python/tests/bend_flux.py +++ b/python/tests/bend_flux.py @@ -33,6 +33,7 @@ def init(self, no_bend=False): mp.Vector3(-0.5 * sx, wvg_ycen + 0.5 * w)] geometry = [mp.Prism(bend_vertices, 0, material=mp.Medium(epsilon=12))] + fcen = 0.15 df = 0.1 sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), component=mp.Ez, @@ -72,32 +73,31 @@ def test_bend_flux(self): fdata = self.sim.get_flux_data(self.refl) expected = [ - (0.1, 3.65231563251e-05, 3.68932495077e-05), - (0.10101010101, 5.55606718876e-05, 5.6065539588e-05), - (0.10202020202, 8.38211697478e-05, 8.44909864736e-05), - (0.10303030303, 0.000125411162229, 0.000126268639045), - (0.10404040404, 0.000186089117531, 0.000187135303398), - (0.105050505051, 0.000273848867869, 0.000275039134667), - (0.106060606061, 0.000399674037745, 0.000400880269423), - (0.107070707071, 0.00057849953593, 0.000579454087881), - (0.108080808081, 0.000830418432986, 0.000830635406881), - (0.109090909091, 0.00118217282661, 0.00118084271347), - (0.110101010101, 0.00166896468348, 0.00166481944189), - (0.111111111111, 0.00233661613864, 0.00232776318321), - (0.112121212121, 0.00324409729096, 0.00322782257917), - (0.113131313131, 0.00446642217385, 0.00443896468822), - (0.114141414141, 0.0060978895019, 0.0060541922825), - (0.115151515152, 0.00825561352398, 0.00818906047274), - (0.116161616162, 0.0110832518495, 0.010985404883), - (0.117171717172, 0.0147547920552, 0.0146151488236), - (0.118181818182, 0.0194782085272, 0.0192840042241), - (0.119191919192, 0.0254987474079, 0.0252348211592), - + (0.09999999999999999, 3.759029112132568e-5, 3.808240613412702e-5), + (0.10101010101010101, 5.718360627275768e-5, 5.787858557476235e-5), + (0.10202020202020202, 8.626178149785403e-5, 8.723095735323989e-5), + (0.10303030303030304, 1.2904302255339786e-4, 1.303734865376113e-4), + (0.10404040404040406, 1.91440977959859e-4, 1.9323229749290062e-4), + (0.10505050505050507, 2.816640015532028e-4, 2.8401853836185545e-4), + (0.10606060606060609, 4.1099104864621264e-4, 4.139936984202774e-4), + (0.10707070707070711, 5.947620615041214e-4, 5.984438388719244e-4), + (0.10808080808080812, 8.536232917474e-4, 8.579019620596675e-4), + (0.10909090909090914, 0.0012150684695463483, 0.0012196587415521445), + (0.11010101010101016, 0.0017153140696826942, 0.0017196000673953989), + (0.11111111111111117, 0.0024015400405837696, 0.0024044040886885415), + (0.11212121212121219, 0.0033345151119150893, 0.0033341192817500924), + (0.11313131313131321, 0.004591607098594056, 0.0045851272360724615), + (0.11414141414141422, 0.006270152592938964, 0.006253468410682328), + (0.11515151515151524, 0.008491126571916682, 0.008458475808916224), + (0.11616161616161626, 0.01140301079799375, 0.011346620857916318), + (0.11717171717171727, 0.015185711272838182, 0.0150954288696048), + (0.11818181818181829, 0.020054321141823886, 0.01991726961986203), + (0.11919191919191931, 0.02626246907430508, 0.026062772997176274), ] res = list(zip(mp.get_flux_freqs(self.trans), mp.get_fluxes(self.trans), mp.get_fluxes(self.refl))) - compare_arrays(self, np.array(expected), np.array(res[:20]), tol=1e-3) + compare_arrays(self, np.array(expected), np.array(res[:20])) # Real run self.sim = None @@ -107,32 +107,30 @@ def test_bend_flux(self): self.sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, self.pt, 1e-3)) expected = [ - (0.09999999999999999, 1.8392235204829767e-5, -7.259467687598002e-6), - (0.10101010101010101, 2.7629932558236724e-5, -1.1107162110079347e-5), - (0.10202020202020202, 4.1001228946782745e-5, -1.687561915798036e-5), - (0.10303030303030304, 6.018966076122556e-5, -2.5425779493709066e-5), - (0.10404040404040406, 8.758554071933231e-5, -3.794958119189475e-5), - (0.10505050505050507, 1.2656696778129198e-4, -5.612512808928115e-5), - (0.10606060606060609, 1.817948859871414e-4, -8.232188174309142e-5), - (0.10707070707070711, 2.594514094902856e-4, -1.1981531280672989e-4), - (0.10808080808080812, 3.6736164837695035e-4, -1.7300125173897737e-4), - (0.10909090909090914, 5.150131339048232e-4, -2.476730940385436e-4), - (0.11010101010101016, 7.136181099374187e-4, -3.5145561406042276e-4), - (0.11111111111111117, 9.76491765781944e-4, -4.944142331545938e-4), - (0.11212121212121219, 0.001320033637882244, -6.897357105189368e-4), - (0.11313131313131321, 0.0017653940714397098, -9.543556354451615e-4), - (0.11414141414141422, 0.0023404727796352857, -0.0013095604571818236), - (0.11515151515151524, 0.0030813962415392098, -0.00178176942635486), - (0.11616161616161626, 0.00403238648982478, -0.0024036650652026112), - (0.11717171717171727, 0.005243320443599316, -0.003215529845495731), - (0.11818181818181829, 0.0067654019326068, -0.004266367104375331), - (0.11919191919191931, 0.008646855439680507, -0.005614491919262783), - + (0.09999999999999999, 2.0170954669155e-5, -7.145323406243266e-6), + (0.10101010101010101, 3.070001955299509e-5, -1.0941385713179106e-5), + (0.10202020202020202, 4.626587530507024e-5, -1.6558923559642594e-5), + (0.10303030303030304, 6.892295567480947e-5, -2.4836579153858833e-5), + (0.10404040404040406, 1.0141364425635142e-4, -3.699982489006669e-5), + (0.10505050505050507, 1.474365085880599e-4, -5.473578454006719e-5), + (0.10606060606060609, 2.1207288913368974e-4, -8.028391176449571e-5), + (0.10707070707070711, 3.0239431640933715e-4, -1.1664290104872828e-4), + (0.10808080808080812, 4.2819888830440315e-4, -1.679368899871326e-4), + (0.10909090909090914, 6.027275359145561e-4, -2.398486110177475e-4), + (0.11010101010101016, 8.431482588216411e-4, -3.39975596024738e-4), + (0.11111111111111117, 0.0011706479532335638, -4.78110538203311e-4), + (0.11212121212121219, 0.0016101926665890486, -6.66666253433329e-4), + (0.11313131313131321, 0.0021903704321079953, -9.214760453657109e-4), + (0.11414141414141422, 0.002944040018564378, -0.0012628973860187313), + (0.11515151515151524, 0.003910493690885223, -0.0017168242626265833), + (0.11616161616161626, 0.005139265971645378, -0.0023153058291112685), + (0.11717171717171727, 0.006694611064215856, -0.0030970004778165016), + (0.11818181818181829, 0.008658482837267306, -0.004108081960927432), + (0.11919191919191931, 0.0111293966280746, -0.005403865420469877), ] - res = list(zip(mp.get_flux_freqs(self.trans), mp.get_fluxes(self.trans), mp.get_fluxes(self.refl))) - compare_arrays(self, np.array(expected), np.array(res[:20]), tol=1e-2) + compare_arrays(self, np.array(expected), np.array(res[:20])) if __name__ == '__main__': diff --git a/scheme/examples/bend-flux.ctl b/scheme/examples/bend-flux.ctl index 9fed645d6..dcaab47dd 100644 --- a/scheme/examples/bend-flux.ctl +++ b/scheme/examples/bend-flux.ctl @@ -16,24 +16,24 @@ (set! geometry (if no-bend? (list - (make block + (make block (center 0 wvg-ycen) - (size infinity w infinity) + (size sx w 0) (material (make dielectric (epsilon 12))))) (list - (make block + (make block (center (* -0.5 pad) wvg-ycen) - (size (- sx pad) w infinity) + (size (- sx pad) w 0) (material (make dielectric (epsilon 12)))) - (make block + (make block (center wvg-xcen (* 0.5 pad)) - (size w (- sy pad) infinity) + (size w (- sy pad) 0) (material (make dielectric (epsilon 12))))))) (define-param fcen 0.15) ; pulse center frequency (define-param df 0.1) ; pulse width (in frequency) (set! sources (list - (make source + (make source (src (make gaussian-src (frequency fcen) (fwidth df))) (component Ez) (center (+ 1 (* -0.5 sx)) wvg-ycen) @@ -52,15 +52,15 @@ (center wvg-xcen (- (/ sy 2) 1.5)) (size (* w 2) 0))))) (define refl ; reflected flux (add-flux fcen df nfreq - (make flux-region + (make flux-region (center (+ (* -0.5 sx) 1.5) wvg-ycen) (size 0 (* w 2))))) ; for normal run, load negated fields to subtract incident from refl. fields (if (not no-bend?) (load-minus-flux "refl-flux" refl)) -(run-sources+ +(run-sources+ (stop-when-fields-decayed 50 Ez - (if no-bend? + (if no-bend? (vector3 (- (/ sx 2) 1.5) wvg-ycen) (vector3 wvg-xcen (- (/ sy 2) 1.5))) 1e-3) From bc4c6a1bb2aeac743f59286c1dded40ef58fdc9e Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Thu, 24 May 2018 07:57:29 -0500 Subject: [PATCH 5/7] Revert "Get new scheme data with fixed size blocks" This reverts commit 0203ab4aa1fc164e9f3b747a84e0b97280bbbd70. --- python/tests/bend_flux.py | 88 ++++++++++++++++++----------------- scheme/examples/bend-flux.ctl | 20 ++++---- 2 files changed, 55 insertions(+), 53 deletions(-) diff --git a/python/tests/bend_flux.py b/python/tests/bend_flux.py index a5272073d..758ce3358 100644 --- a/python/tests/bend_flux.py +++ b/python/tests/bend_flux.py @@ -33,7 +33,6 @@ def init(self, no_bend=False): mp.Vector3(-0.5 * sx, wvg_ycen + 0.5 * w)] geometry = [mp.Prism(bend_vertices, 0, material=mp.Medium(epsilon=12))] - fcen = 0.15 df = 0.1 sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), component=mp.Ez, @@ -73,31 +72,32 @@ def test_bend_flux(self): fdata = self.sim.get_flux_data(self.refl) expected = [ - (0.09999999999999999, 3.759029112132568e-5, 3.808240613412702e-5), - (0.10101010101010101, 5.718360627275768e-5, 5.787858557476235e-5), - (0.10202020202020202, 8.626178149785403e-5, 8.723095735323989e-5), - (0.10303030303030304, 1.2904302255339786e-4, 1.303734865376113e-4), - (0.10404040404040406, 1.91440977959859e-4, 1.9323229749290062e-4), - (0.10505050505050507, 2.816640015532028e-4, 2.8401853836185545e-4), - (0.10606060606060609, 4.1099104864621264e-4, 4.139936984202774e-4), - (0.10707070707070711, 5.947620615041214e-4, 5.984438388719244e-4), - (0.10808080808080812, 8.536232917474e-4, 8.579019620596675e-4), - (0.10909090909090914, 0.0012150684695463483, 0.0012196587415521445), - (0.11010101010101016, 0.0017153140696826942, 0.0017196000673953989), - (0.11111111111111117, 0.0024015400405837696, 0.0024044040886885415), - (0.11212121212121219, 0.0033345151119150893, 0.0033341192817500924), - (0.11313131313131321, 0.004591607098594056, 0.0045851272360724615), - (0.11414141414141422, 0.006270152592938964, 0.006253468410682328), - (0.11515151515151524, 0.008491126571916682, 0.008458475808916224), - (0.11616161616161626, 0.01140301079799375, 0.011346620857916318), - (0.11717171717171727, 0.015185711272838182, 0.0150954288696048), - (0.11818181818181829, 0.020054321141823886, 0.01991726961986203), - (0.11919191919191931, 0.02626246907430508, 0.026062772997176274), + (0.1, 3.65231563251e-05, 3.68932495077e-05), + (0.10101010101, 5.55606718876e-05, 5.6065539588e-05), + (0.10202020202, 8.38211697478e-05, 8.44909864736e-05), + (0.10303030303, 0.000125411162229, 0.000126268639045), + (0.10404040404, 0.000186089117531, 0.000187135303398), + (0.105050505051, 0.000273848867869, 0.000275039134667), + (0.106060606061, 0.000399674037745, 0.000400880269423), + (0.107070707071, 0.00057849953593, 0.000579454087881), + (0.108080808081, 0.000830418432986, 0.000830635406881), + (0.109090909091, 0.00118217282661, 0.00118084271347), + (0.110101010101, 0.00166896468348, 0.00166481944189), + (0.111111111111, 0.00233661613864, 0.00232776318321), + (0.112121212121, 0.00324409729096, 0.00322782257917), + (0.113131313131, 0.00446642217385, 0.00443896468822), + (0.114141414141, 0.0060978895019, 0.0060541922825), + (0.115151515152, 0.00825561352398, 0.00818906047274), + (0.116161616162, 0.0110832518495, 0.010985404883), + (0.117171717172, 0.0147547920552, 0.0146151488236), + (0.118181818182, 0.0194782085272, 0.0192840042241), + (0.119191919192, 0.0254987474079, 0.0252348211592), + ] res = list(zip(mp.get_flux_freqs(self.trans), mp.get_fluxes(self.trans), mp.get_fluxes(self.refl))) - compare_arrays(self, np.array(expected), np.array(res[:20])) + compare_arrays(self, np.array(expected), np.array(res[:20]), tol=1e-3) # Real run self.sim = None @@ -107,30 +107,32 @@ def test_bend_flux(self): self.sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, self.pt, 1e-3)) expected = [ - (0.09999999999999999, 2.0170954669155e-5, -7.145323406243266e-6), - (0.10101010101010101, 3.070001955299509e-5, -1.0941385713179106e-5), - (0.10202020202020202, 4.626587530507024e-5, -1.6558923559642594e-5), - (0.10303030303030304, 6.892295567480947e-5, -2.4836579153858833e-5), - (0.10404040404040406, 1.0141364425635142e-4, -3.699982489006669e-5), - (0.10505050505050507, 1.474365085880599e-4, -5.473578454006719e-5), - (0.10606060606060609, 2.1207288913368974e-4, -8.028391176449571e-5), - (0.10707070707070711, 3.0239431640933715e-4, -1.1664290104872828e-4), - (0.10808080808080812, 4.2819888830440315e-4, -1.679368899871326e-4), - (0.10909090909090914, 6.027275359145561e-4, -2.398486110177475e-4), - (0.11010101010101016, 8.431482588216411e-4, -3.39975596024738e-4), - (0.11111111111111117, 0.0011706479532335638, -4.78110538203311e-4), - (0.11212121212121219, 0.0016101926665890486, -6.66666253433329e-4), - (0.11313131313131321, 0.0021903704321079953, -9.214760453657109e-4), - (0.11414141414141422, 0.002944040018564378, -0.0012628973860187313), - (0.11515151515151524, 0.003910493690885223, -0.0017168242626265833), - (0.11616161616161626, 0.005139265971645378, -0.0023153058291112685), - (0.11717171717171727, 0.006694611064215856, -0.0030970004778165016), - (0.11818181818181829, 0.008658482837267306, -0.004108081960927432), - (0.11919191919191931, 0.0111293966280746, -0.005403865420469877), + (0.09999999999999999, 1.8392235204829767e-5, -7.259467687598002e-6), + (0.10101010101010101, 2.7629932558236724e-5, -1.1107162110079347e-5), + (0.10202020202020202, 4.1001228946782745e-5, -1.687561915798036e-5), + (0.10303030303030304, 6.018966076122556e-5, -2.5425779493709066e-5), + (0.10404040404040406, 8.758554071933231e-5, -3.794958119189475e-5), + (0.10505050505050507, 1.2656696778129198e-4, -5.612512808928115e-5), + (0.10606060606060609, 1.817948859871414e-4, -8.232188174309142e-5), + (0.10707070707070711, 2.594514094902856e-4, -1.1981531280672989e-4), + (0.10808080808080812, 3.6736164837695035e-4, -1.7300125173897737e-4), + (0.10909090909090914, 5.150131339048232e-4, -2.476730940385436e-4), + (0.11010101010101016, 7.136181099374187e-4, -3.5145561406042276e-4), + (0.11111111111111117, 9.76491765781944e-4, -4.944142331545938e-4), + (0.11212121212121219, 0.001320033637882244, -6.897357105189368e-4), + (0.11313131313131321, 0.0017653940714397098, -9.543556354451615e-4), + (0.11414141414141422, 0.0023404727796352857, -0.0013095604571818236), + (0.11515151515151524, 0.0030813962415392098, -0.00178176942635486), + (0.11616161616161626, 0.00403238648982478, -0.0024036650652026112), + (0.11717171717171727, 0.005243320443599316, -0.003215529845495731), + (0.11818181818181829, 0.0067654019326068, -0.004266367104375331), + (0.11919191919191931, 0.008646855439680507, -0.005614491919262783), + ] + res = list(zip(mp.get_flux_freqs(self.trans), mp.get_fluxes(self.trans), mp.get_fluxes(self.refl))) - compare_arrays(self, np.array(expected), np.array(res[:20])) + compare_arrays(self, np.array(expected), np.array(res[:20]), tol=1e-2) if __name__ == '__main__': diff --git a/scheme/examples/bend-flux.ctl b/scheme/examples/bend-flux.ctl index dcaab47dd..9fed645d6 100644 --- a/scheme/examples/bend-flux.ctl +++ b/scheme/examples/bend-flux.ctl @@ -16,24 +16,24 @@ (set! geometry (if no-bend? (list - (make block + (make block (center 0 wvg-ycen) - (size sx w 0) + (size infinity w infinity) (material (make dielectric (epsilon 12))))) (list - (make block + (make block (center (* -0.5 pad) wvg-ycen) - (size (- sx pad) w 0) + (size (- sx pad) w infinity) (material (make dielectric (epsilon 12)))) - (make block + (make block (center wvg-xcen (* 0.5 pad)) - (size w (- sy pad) 0) + (size w (- sy pad) infinity) (material (make dielectric (epsilon 12))))))) (define-param fcen 0.15) ; pulse center frequency (define-param df 0.1) ; pulse width (in frequency) (set! sources (list - (make source + (make source (src (make gaussian-src (frequency fcen) (fwidth df))) (component Ez) (center (+ 1 (* -0.5 sx)) wvg-ycen) @@ -52,15 +52,15 @@ (center wvg-xcen (- (/ sy 2) 1.5)) (size (* w 2) 0))))) (define refl ; reflected flux (add-flux fcen df nfreq - (make flux-region + (make flux-region (center (+ (* -0.5 sx) 1.5) wvg-ycen) (size 0 (* w 2))))) ; for normal run, load negated fields to subtract incident from refl. fields (if (not no-bend?) (load-minus-flux "refl-flux" refl)) -(run-sources+ +(run-sources+ (stop-when-fields-decayed 50 Ez - (if no-bend? + (if no-bend? (vector3 (- (/ sx 2) 1.5) wvg-ycen) (vector3 wvg-xcen (- (/ sy 2) 1.5))) 1e-3) From efeeeac5a5f8c59b89ae82eb24627eb36c2561e2 Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Thu, 24 May 2018 08:18:02 -0500 Subject: [PATCH 6/7] height is double --- python/meep.i | 13 +++++++++++++ python/typemap_utils.cpp | 17 ++--------------- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/python/meep.i b/python/meep.i index 929d56a9b..5ee7ff4ed 100644 --- a/python/meep.i +++ b/python/meep.i @@ -58,6 +58,19 @@ typedef struct { #include "typemap_utils.cpp" +static int get_attr_int(PyObject *py_obj, int *result, const char *name) { + PyObject *py_attr = PyObject_GetAttrString(py_obj, name); + + if (!py_attr) { + PyErr_Format(PyExc_ValueError, "Class attribute '%s' is None\n", name); + return 0; + } + + *result = PyInteger_AsLong(py_attr); + Py_XDECREF(py_attr); + return 1; +} + static PyObject *py_source_time_object() { static PyObject *source_time_object = NULL; if (source_time_object == NULL) { diff --git a/python/typemap_utils.cpp b/python/typemap_utils.cpp index fd695fc43..69f1e94ef 100644 --- a/python/typemap_utils.cpp +++ b/python/typemap_utils.cpp @@ -263,19 +263,6 @@ static int get_attr_dbl(PyObject *py_obj, double *result, const char *name) { return 1; } -static int get_attr_int(PyObject *py_obj, int *result, const char *name) { - PyObject *py_attr = PyObject_GetAttrString(py_obj, name); - - if (!py_attr) { - PyErr_Format(PyExc_ValueError, "Class attribute '%s' is None\n", name); - return 0; - } - - *result = PyInteger_AsLong(py_attr); - Py_XDECREF(py_attr); - return 1; -} - static int get_attr_material(PyObject *po, material_type *m) { PyObject *py_material = PyObject_GetAttrString(po, "material"); @@ -556,11 +543,11 @@ static int pyellipsoid_to_ellipsoid(PyObject *py_ell, geometric_object *e) { static int pyprism_to_prism(PyObject *py_prism, geometric_object *p) { material_type material; - int height; + double height; vector3 axis; if (!get_attr_material(py_prism, &material) || - !get_attr_int(py_prism, &height, "height") || + !get_attr_dbl(py_prism, &height, "height") || !get_attr_v3(py_prism, &axis, "axis")) { return 0; From 000d8c337645fe6d1bb6e766ab109e668bb7ad9b Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Fri, 1 Jun 2018 08:45:41 -0500 Subject: [PATCH 7/7] Extend prism out of cell and increase height. --- python/tests/bend_flux.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/python/tests/bend_flux.py b/python/tests/bend_flux.py index 758ce3358..f351edf46 100644 --- a/python/tests/bend_flux.py +++ b/python/tests/bend_flux.py @@ -16,14 +16,15 @@ def init(self, no_bend=False): w = 1 wvg_ycen = -0.5 * (sy - w - (2 * pad)) wvg_xcen = 0.5 * (sx - w - (2 * pad)) + height = 100 if no_bend: - no_bend_vertices = [mp.Vector3(-0.5 * sx, wvg_ycen - 0.5 * w), - mp.Vector3(+0.5 * sx, wvg_ycen - 0.5 * w), - mp.Vector3(+0.5 * sx, wvg_ycen + 0.5 * w), - mp.Vector3(-0.5 * sx, wvg_ycen + 0.5 * w)] + no_bend_vertices = [mp.Vector3(-0.5 * sx - 5, wvg_ycen - 0.5 * w), + mp.Vector3(+0.5 * sx + 5, wvg_ycen - 0.5 * w), + mp.Vector3(+0.5 * sx + 5, wvg_ycen + 0.5 * w), + mp.Vector3(-0.5 * sx - 5, wvg_ycen + 0.5 * w)] - geometry = [mp.Prism(no_bend_vertices, 0, material=mp.Medium(epsilon=12))] + geometry = [mp.Prism(no_bend_vertices, height, material=mp.Medium(epsilon=12))] else: bend_vertices = [mp.Vector3(-0.5 * sx, wvg_ycen - 0.5 * w), mp.Vector3(wvg_xcen + 0.5 * w, wvg_ycen - 0.5 * w), @@ -32,7 +33,8 @@ def init(self, no_bend=False): mp.Vector3(wvg_xcen - 0.5 * w, wvg_ycen + 0.5 * w), mp.Vector3(-0.5 * sx, wvg_ycen + 0.5 * w)] - geometry = [mp.Prism(bend_vertices, 0, material=mp.Medium(epsilon=12))] + geometry = [mp.Prism(bend_vertices, height, material=mp.Medium(epsilon=12))] + fcen = 0.15 df = 0.1 sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), component=mp.Ez, @@ -92,12 +94,11 @@ def test_bend_flux(self): (0.117171717172, 0.0147547920552, 0.0146151488236), (0.118181818182, 0.0194782085272, 0.0192840042241), (0.119191919192, 0.0254987474079, 0.0252348211592), - ] res = list(zip(mp.get_flux_freqs(self.trans), mp.get_fluxes(self.trans), mp.get_fluxes(self.refl))) - compare_arrays(self, np.array(expected), np.array(res[:20]), tol=1e-3) + compare_arrays(self, np.array(expected), np.array(res[:20]), tol=1e-7) # Real run self.sim = None @@ -132,7 +133,7 @@ def test_bend_flux(self): res = list(zip(mp.get_flux_freqs(self.trans), mp.get_fluxes(self.trans), mp.get_fluxes(self.refl))) - compare_arrays(self, np.array(expected), np.array(res[:20]), tol=1e-2) + compare_arrays(self, np.array(expected), np.array(res[:20]), tol=1e-3) if __name__ == '__main__':