From 83fd9ef41c071b57783ff58fb0a5d6eda7bceb29 Mon Sep 17 00:00:00 2001 From: smartalecH Date: Fri, 12 Aug 2022 10:51:10 -0400 Subject: [PATCH 1/5] add function to update m number --- python/adjoint/optimization_problem.py | 6 ++-- python/simulation.py | 19 ++++++++++++ src/meepgeom.cpp | 40 +++++++++++++++----------- 3 files changed, 45 insertions(+), 20 deletions(-) diff --git a/python/adjoint/optimization_problem.py b/python/adjoint/optimization_problem.py index 707eb8beb..a7899580e 100644 --- a/python/adjoint/optimization_problem.py +++ b/python/adjoint/optimization_problem.py @@ -226,7 +226,7 @@ def adjoint_run(self): # flip the m number if utils._check_if_cylindrical(self.sim): - self.sim.m = -self.sim.m + self.sim.change_m_number(-self.sim.m) # flip the k point if self.sim.k_point: @@ -261,11 +261,11 @@ def adjoint_run(self): # reset the m number if utils._check_if_cylindrical(self.sim): - self.sim.m = -self.sim.m + self.sim.change_m_number(-self.sim.m) # reset the k point if self.sim.k_point: - self.sim.k_point *= -1 + self.sim.change_k_point(-1 * self.sim.k_point) # update optimizer's state self.current_state = "ADJ" diff --git a/python/simulation.py b/python/simulation.py index a8afd9e3a..9662d4f5b 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -4278,6 +4278,25 @@ def change_k_point(self, k): self.fields.use_bloch( py_v3_to_vec(self.dimensions, self.k_point, self.is_cylindrical) ) + + def change_m_number(self, m): + """ + Change the `m` number (the angular ϕ dependence). + """ + self.m = m + + if self.fields: + needs_complex_fields = not ( + not self.m or self.m == 0 + ) + + if needs_complex_fields and self.fields.is_real: + self.fields = None + self._is_initialized = False + self.init_sim() + else: + if self.m is not None: + self.fields.m = self.m def change_sources(self, new_sources): """ diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp index 6b83d109b..a37ab2442 100644 --- a/src/meepgeom.cpp +++ b/src/meepgeom.cpp @@ -2651,16 +2651,25 @@ get_material_gradient(const meep::vec &r, // current point material_type md; geps->get_material_pt(md, r); + // get the tensor column index corresponding to the forward component int dir_idx = 0; - if (forward_c == meep::Dx || forward_c == meep::Dr) - dir_idx = 0; - else if (forward_c == meep::Dy || forward_c == meep::Dp) - dir_idx = 1; - else if (forward_c == meep::Dz) - dir_idx = 2; - else - meep::abort("Invalid adjoint field component"); + switch (meep::component_direction(forward_c)) { + case meep::X: + case meep::R: + dir_idx = 0; + break; + case meep::Y: + case meep::P: + dir_idx = 1; + break; + case meep::Z: + dir_idx = 2; + break; + case meep::NO_DIRECTION: + meep::abort("Invalid forward component!\n"); + } + // materials are non-dispersive if (md->trivial) { const double sd = 1.0; // FIXME: make user-changable? meep::volume v(r); @@ -2668,18 +2677,17 @@ get_material_gradient(const meep::vec &r, // current point v.set_direction_min(d, r.in_direction(d) - 0.5 * gv.inva * sd); v.set_direction_max(d, r.in_direction(d) + 0.5 * gv.inva * sd); } - double row_1[3], row_2[3], dA_du[3]; + double row_1[3], row_2[3]; double orig = u[idx]; u[idx] -= du; geps->eff_chi1inv_row(adjoint_c, row_1, v, geps->tol, geps->maxeval); u[idx] += 2 * du; geps->eff_chi1inv_row(adjoint_c, row_2, v, geps->tol, geps->maxeval); u[idx] = orig; + return fields_f * (row_1[dir_idx] - row_2[dir_idx]) / (2 * du); - for (int i = 0; i < 3; i++) - dA_du[i] = (row_1[i] - row_2[i]) / (2 * du); - return dA_du[dir_idx] * fields_f; } + // materials have some dispersion else { double orig = u[idx]; std::complex row_1[3], row_2[3], dA_du[3]; @@ -2689,9 +2697,7 @@ get_material_gradient(const meep::vec &r, // current point eff_chi1inv_row_disp(adjoint_c, row_2, r, freq, geps); u[idx] = orig; - for (int i = 0; i < 3; i++) - dA_du[i] = (row_1[i] - row_2[i]) / (2 * du); - return dA_du[dir_idx] * fields_f * cond_cmp(forward_c, r, freq, geps); + return fields_f * (row_1[dir_idx] - row_2[dir_idx]) / (2 * du) * cond_cmp(forward_c, r, freq, geps); } } @@ -2720,8 +2726,8 @@ void add_interpolate_weights(double rx, double ry, double rz, double *data, int /* define a macro to give us data(x,y,z) on the grid, in row-major order (the order used by HDF5): */ #define IDX(x, y, z) (((x)*ny + (y)) * nz + (z)) * stride -#define D(x, y, z) (data[(((x)*ny + (y)) * nz + (z)) * stride]) -#define U(x, y, z) (udata[(((x)*ny + (y)) * nz + (z)) * stride]) +#define D(x, y, z) (data[IDX(x, y, z)]) +#define U(x, y, z) (udata[IDX(x, y, z)]) u = (((U(x1, y1, z1) * (1.0 - dx) + U(x2, y1, z1) * dx) * (1.0 - dy) + (U(x1, y2, z1) * (1.0 - dx) + U(x2, y2, z1) * dx) * dy) * From 2dff39ef02a81431d36e566725aa8e04f12327fc Mon Sep 17 00:00:00 2001 From: smartalecH Date: Fri, 12 Aug 2022 13:58:16 -0400 Subject: [PATCH 2/5] update change m number too --- python/simulation.py | 4 ++-- src/fields.cpp | 11 +++++++++++ src/meep.hpp | 2 ++ 3 files changed, 15 insertions(+), 2 deletions(-) diff --git a/python/simulation.py b/python/simulation.py index 9662d4f5b..13aecfbe0 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -4284,7 +4284,7 @@ def change_m_number(self, m): Change the `m` number (the angular ϕ dependence). """ self.m = m - + if self.fields: needs_complex_fields = not ( not self.m or self.m == 0 @@ -4296,7 +4296,7 @@ def change_m_number(self, m): self.init_sim() else: if self.m is not None: - self.fields.m = self.m + self.fields.change_m_number(m) def change_sources(self, new_sources): """ diff --git a/src/fields.cpp b/src/fields.cpp index b9f43f338..a8c568fd3 100644 --- a/src/fields.cpp +++ b/src/fields.cpp @@ -793,4 +793,15 @@ bool operator==(const comms_key &lhs, const comms_key &rhs) { return (lhs.ft == rhs.ft) && (lhs.phase == rhs.phase) && (lhs.pair == rhs.pair); } +void fields::change_m_number(double new_m){ + m = new_m; + for (int i = 0; i < num_chunks; i++) { + chunks[i]->change_m_number(new_m); + } +} + +void fields_chunk::change_m_number(double new_m) { + m = new_m; +} + } // namespace meep diff --git a/src/meep.hpp b/src/meep.hpp index f5a721e81..9b82072a2 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1546,6 +1546,7 @@ class fields_chunk { void remove_sources(); void remove_susceptibilities(bool shared_chunks); void zero_fields(); + void change_m_number(double new_m); // update_eh.cpp bool needs_W_prev(component c) const; @@ -1751,6 +1752,7 @@ class fields { void remove_fluxes(); void reset(); void log(const char *prefix = ""); + void change_m_number(double new_m); // time.cpp std::vector time_spent_on(time_sink sink); From 2230adbf1bc5a7d2cf6b11ce03549cbde1a9bd5f Mon Sep 17 00:00:00 2001 From: smartalecH Date: Fri, 12 Aug 2022 14:48:26 -0400 Subject: [PATCH 3/5] precommit cleanup --- python/simulation.py | 6 ++---- src/fields.cpp | 6 ++---- src/meepgeom.cpp | 23 ++++++++--------------- 3 files changed, 12 insertions(+), 23 deletions(-) diff --git a/python/simulation.py b/python/simulation.py index 13aecfbe0..6c3e15d63 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -4278,7 +4278,7 @@ def change_k_point(self, k): self.fields.use_bloch( py_v3_to_vec(self.dimensions, self.k_point, self.is_cylindrical) ) - + def change_m_number(self, m): """ Change the `m` number (the angular ϕ dependence). @@ -4286,9 +4286,7 @@ def change_m_number(self, m): self.m = m if self.fields: - needs_complex_fields = not ( - not self.m or self.m == 0 - ) + needs_complex_fields = not (not self.m or self.m == 0) if needs_complex_fields and self.fields.is_real: self.fields = None diff --git a/src/fields.cpp b/src/fields.cpp index a8c568fd3..6d16eb748 100644 --- a/src/fields.cpp +++ b/src/fields.cpp @@ -793,15 +793,13 @@ bool operator==(const comms_key &lhs, const comms_key &rhs) { return (lhs.ft == rhs.ft) && (lhs.phase == rhs.phase) && (lhs.pair == rhs.pair); } -void fields::change_m_number(double new_m){ +void fields::change_m_number(double new_m) { m = new_m; for (int i = 0; i < num_chunks; i++) { chunks[i]->change_m_number(new_m); } } -void fields_chunk::change_m_number(double new_m) { - m = new_m; -} +void fields_chunk::change_m_number(double new_m) { m = new_m; } } // namespace meep diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp index a37ab2442..15fa04076 100644 --- a/src/meepgeom.cpp +++ b/src/meepgeom.cpp @@ -2654,19 +2654,12 @@ get_material_gradient(const meep::vec &r, // current point // get the tensor column index corresponding to the forward component int dir_idx = 0; switch (meep::component_direction(forward_c)) { - case meep::X: - case meep::R: - dir_idx = 0; - break; - case meep::Y: - case meep::P: - dir_idx = 1; - break; - case meep::Z: - dir_idx = 2; - break; - case meep::NO_DIRECTION: - meep::abort("Invalid forward component!\n"); + case meep::X: + case meep::R: dir_idx = 0; break; + case meep::Y: + case meep::P: dir_idx = 1; break; + case meep::Z: dir_idx = 2; break; + case meep::NO_DIRECTION: meep::abort("Invalid forward component!\n"); } // materials are non-dispersive @@ -2685,7 +2678,6 @@ get_material_gradient(const meep::vec &r, // current point geps->eff_chi1inv_row(adjoint_c, row_2, v, geps->tol, geps->maxeval); u[idx] = orig; return fields_f * (row_1[dir_idx] - row_2[dir_idx]) / (2 * du); - } // materials have some dispersion else { @@ -2697,7 +2689,8 @@ get_material_gradient(const meep::vec &r, // current point eff_chi1inv_row_disp(adjoint_c, row_2, r, freq, geps); u[idx] = orig; - return fields_f * (row_1[dir_idx] - row_2[dir_idx]) / (2 * du) * cond_cmp(forward_c, r, freq, geps); + return fields_f * (row_1[dir_idx] - row_2[dir_idx]) / (2 * du) * + cond_cmp(forward_c, r, freq, geps); } } From 514128e1118a427206285311861a647e97584c90 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Thu, 22 Sep 2022 09:56:24 -0700 Subject: [PATCH 4/5] naming fixes, added check to change_m(), improved tests --- python/adjoint/optimization_problem.py | 4 ++-- python/simulation.py | 7 +++---- python/tests/test_adjoint_cyl.py | 24 +++++++++++++++--------- src/fields.cpp | 14 +++++++++++--- src/meep.hpp | 4 ++-- 5 files changed, 33 insertions(+), 20 deletions(-) diff --git a/python/adjoint/optimization_problem.py b/python/adjoint/optimization_problem.py index a7899580e..8e137144b 100644 --- a/python/adjoint/optimization_problem.py +++ b/python/adjoint/optimization_problem.py @@ -226,7 +226,7 @@ def adjoint_run(self): # flip the m number if utils._check_if_cylindrical(self.sim): - self.sim.change_m_number(-self.sim.m) + self.sim.change_m(-self.sim.m) # flip the k point if self.sim.k_point: @@ -261,7 +261,7 @@ def adjoint_run(self): # reset the m number if utils._check_if_cylindrical(self.sim): - self.sim.change_m_number(-self.sim.m) + self.sim.change_m(-self.sim.m) # reset the k point if self.sim.k_point: diff --git a/python/simulation.py b/python/simulation.py index 6c3e15d63..d0c5ca407 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -4279,9 +4279,8 @@ def change_k_point(self, k): py_v3_to_vec(self.dimensions, self.k_point, self.is_cylindrical) ) - def change_m_number(self, m): - """ - Change the `m` number (the angular ϕ dependence). + def change_m(self, m: float)->None: + """Changes the simulation's `m` number (the angular ϕ dependence). """ self.m = m @@ -4294,7 +4293,7 @@ def change_m_number(self, m): self.init_sim() else: if self.m is not None: - self.fields.change_m_number(m) + self.fields.change_m(m) def change_sources(self, new_sources): """ diff --git a/python/tests/test_adjoint_cyl.py b/python/tests/test_adjoint_cyl.py index 98bc01d33..69b8e64c9 100644 --- a/python/tests/test_adjoint_cyl.py +++ b/python/tests/test_adjoint_cyl.py @@ -12,11 +12,11 @@ from autograd import numpy as npa from autograd import tensor_jacobian_product from utils import ApproxComparisonTestCase +import parameterized rng = np.random.RandomState(2) resolution = 20 dimensions = mp.CYLINDRICAL -m = 0 Si = mp.Medium(index=3.4) SiO2 = mp.Medium(index=1.44) @@ -48,7 +48,7 @@ dp = deps * rng.rand(Nr * Nz) -def forward_simulation(design_params): +def forward_simulation(design_params,m,far_x): matgrid = mp.MaterialGrid( mp.Vector3(Nr, 0, Nz), SiO2, Si, weights=design_params.reshape(Nr, 1, Nz) ) @@ -72,7 +72,7 @@ def forward_simulation(design_params): ) frequencies = [fcen] - far_x = [mp.Vector3(5, 0, 20)] + mode = sim.add_near2far( frequencies, mp.Near2FarRegion( @@ -89,7 +89,7 @@ def forward_simulation(design_params): return abs(Er[0]) ** 2 -def adjoint_solver(design_params): +def adjoint_solver(design_params,m,far_x): design_variables = mp.MaterialGrid(mp.Vector3(Nr, 0, Nz), SiO2, Si) design_region = mpa.DesignRegion( @@ -117,7 +117,6 @@ def adjoint_solver(design_params): m=m, ) - far_x = [mp.Vector3(5, 0, 20)] NearRegions = [ mp.Near2FarRegion( center=mp.Vector3(design_r / 2, 0, (sz / 2 - dpml + design_z / 2) / 2), @@ -149,12 +148,19 @@ def J(alpha): class TestAdjointSolver(ApproxComparisonTestCase): - def test_adjoint_solver_cyl_n2f_fields(self): + @parameterized.parameterized.expand([ + (0,[mp.Vector3(5, 0, 20)]), + (0,[mp.Vector3(4, 0, 28)]), + (-1,[mp.Vector3(5, 0, 20)]), + (1.2,[mp.Vector3(5, 0, 20)]) + ]) + def test_adjoint_solver_cyl_n2f_fields(self,m,far_x): print("*** TESTING CYLINDRICAL Near2Far ADJOINT FEATURES ***") - adjsol_obj, adjsol_grad = adjoint_solver(p) + print("Current test: m={}, far_x={}".format(m,far_x)) + adjsol_obj, adjsol_grad = adjoint_solver(p,m,far_x) ## compute unperturbed S12 - S12_unperturbed = forward_simulation(p) + S12_unperturbed = forward_simulation(p,m,far_x) ## compare objective results print( @@ -164,7 +170,7 @@ def test_adjoint_solver_cyl_n2f_fields(self): self.assertClose(adjsol_obj, S12_unperturbed, epsilon=1e-3) ## compute perturbed S12 - S12_perturbed = forward_simulation(p + dp) + S12_perturbed = forward_simulation(p + dp,m,far_x) ## compare gradients if adjsol_grad.ndim < 2: diff --git a/src/fields.cpp b/src/fields.cpp index 6d16eb748..59d1299dc 100644 --- a/src/fields.cpp +++ b/src/fields.cpp @@ -793,13 +793,21 @@ bool operator==(const comms_key &lhs, const comms_key &rhs) { return (lhs.ft == rhs.ft) && (lhs.phase == rhs.phase) && (lhs.pair == rhs.pair); } -void fields::change_m_number(double new_m) { +void fields::change_m(double new_m) { m = new_m; + if ((new_m != 0) && (is_real)){ + meep::abort("The simulation must be reinitialized if switching to complex fields!\n"); + } + + if ((new_m == 0) && (!is_real)){ + use_real_fields(); + } + for (int i = 0; i < num_chunks; i++) { - chunks[i]->change_m_number(new_m); + chunks[i]->change_m(new_m); } } -void fields_chunk::change_m_number(double new_m) { m = new_m; } +void fields_chunk::change_m(double new_m) { m = new_m; } } // namespace meep diff --git a/src/meep.hpp b/src/meep.hpp index 9b82072a2..d7e9a2ce3 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1546,7 +1546,7 @@ class fields_chunk { void remove_sources(); void remove_susceptibilities(bool shared_chunks); void zero_fields(); - void change_m_number(double new_m); + void change_m(double new_m); // update_eh.cpp bool needs_W_prev(component c) const; @@ -1752,7 +1752,7 @@ class fields { void remove_fluxes(); void reset(); void log(const char *prefix = ""); - void change_m_number(double new_m); + void change_m(double new_m); // time.cpp std::vector time_spent_on(time_sink sink); From ece632981528da8a3c63425292937c93db56ec8a Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Thu, 22 Sep 2022 10:00:22 -0700 Subject: [PATCH 5/5] fix precommit errors --- python/simulation.py | 5 ++--- python/tests/test_adjoint_cyl.py | 28 +++++++++++++++------------- src/fields.cpp | 6 ++---- 3 files changed, 19 insertions(+), 20 deletions(-) diff --git a/python/simulation.py b/python/simulation.py index d0c5ca407..fe7c38ed3 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -4279,9 +4279,8 @@ def change_k_point(self, k): py_v3_to_vec(self.dimensions, self.k_point, self.is_cylindrical) ) - def change_m(self, m: float)->None: - """Changes the simulation's `m` number (the angular ϕ dependence). - """ + def change_m(self, m: float) -> None: + """Changes the simulation's `m` number (the angular ϕ dependence).""" self.m = m if self.fields: diff --git a/python/tests/test_adjoint_cyl.py b/python/tests/test_adjoint_cyl.py index 69b8e64c9..2ff6c0884 100644 --- a/python/tests/test_adjoint_cyl.py +++ b/python/tests/test_adjoint_cyl.py @@ -48,7 +48,7 @@ dp = deps * rng.rand(Nr * Nz) -def forward_simulation(design_params,m,far_x): +def forward_simulation(design_params, m, far_x): matgrid = mp.MaterialGrid( mp.Vector3(Nr, 0, Nz), SiO2, Si, weights=design_params.reshape(Nr, 1, Nz) ) @@ -89,7 +89,7 @@ def forward_simulation(design_params,m,far_x): return abs(Er[0]) ** 2 -def adjoint_solver(design_params,m,far_x): +def adjoint_solver(design_params, m, far_x): design_variables = mp.MaterialGrid(mp.Vector3(Nr, 0, Nz), SiO2, Si) design_region = mpa.DesignRegion( @@ -148,19 +148,21 @@ def J(alpha): class TestAdjointSolver(ApproxComparisonTestCase): - @parameterized.parameterized.expand([ - (0,[mp.Vector3(5, 0, 20)]), - (0,[mp.Vector3(4, 0, 28)]), - (-1,[mp.Vector3(5, 0, 20)]), - (1.2,[mp.Vector3(5, 0, 20)]) - ]) - def test_adjoint_solver_cyl_n2f_fields(self,m,far_x): + @parameterized.parameterized.expand( + [ + (0, [mp.Vector3(5, 0, 20)]), + (0, [mp.Vector3(4, 0, 28)]), + (-1, [mp.Vector3(5, 0, 20)]), + (1.2, [mp.Vector3(5, 0, 20)]), + ] + ) + def test_adjoint_solver_cyl_n2f_fields(self, m, far_x): print("*** TESTING CYLINDRICAL Near2Far ADJOINT FEATURES ***") - print("Current test: m={}, far_x={}".format(m,far_x)) - adjsol_obj, adjsol_grad = adjoint_solver(p,m,far_x) + print(f"Current test: m={m}, far_x={far_x}") + adjsol_obj, adjsol_grad = adjoint_solver(p, m, far_x) ## compute unperturbed S12 - S12_unperturbed = forward_simulation(p,m,far_x) + S12_unperturbed = forward_simulation(p, m, far_x) ## compare objective results print( @@ -170,7 +172,7 @@ def test_adjoint_solver_cyl_n2f_fields(self,m,far_x): self.assertClose(adjsol_obj, S12_unperturbed, epsilon=1e-3) ## compute perturbed S12 - S12_perturbed = forward_simulation(p + dp,m,far_x) + S12_perturbed = forward_simulation(p + dp, m, far_x) ## compare gradients if adjsol_grad.ndim < 2: diff --git a/src/fields.cpp b/src/fields.cpp index 59d1299dc..2eee98ea8 100644 --- a/src/fields.cpp +++ b/src/fields.cpp @@ -795,13 +795,11 @@ bool operator==(const comms_key &lhs, const comms_key &rhs) { void fields::change_m(double new_m) { m = new_m; - if ((new_m != 0) && (is_real)){ + if ((new_m != 0) && (is_real)) { meep::abort("The simulation must be reinitialized if switching to complex fields!\n"); } - if ((new_m == 0) && (!is_real)){ - use_real_fields(); - } + if ((new_m == 0) && (!is_real)) { use_real_fields(); } for (int i = 0; i < num_chunks; i++) { chunks[i]->change_m(new_m);