Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Resolve m-number adjoint bugs introduced by 1855 #2194

Merged
merged 5 commits into from
Sep 22, 2022
Merged

Resolve m-number adjoint bugs introduced by 1855 #2194

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
83fd9ef
add function to update m number
smartalecH Aug 12, 2022
2dff39e
update change m number too
smartalecH Aug 12, 2022
2230adb
precommit cleanup
smartalecH Aug 12, 2022
514128e
naming fixes, added check to change_m(), improved tests
Sep 22, 2022
ece6329
fix precommit errors
Sep 22, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions python/adjoint/optimization_problem.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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"
Expand Down
17 changes: 17 additions & 0 deletions python/simulation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4279,6 +4279,23 @@ def change_k_point(self, k):
py_v3_to_vec(self.dimensions, self.k_point, self.is_cylindrical)
)

def change_m_number(self, m):
smartalecH marked this conversation as resolved.
Show resolved Hide resolved
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe just call it change_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()
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would be better to just add a boolean argument to fields::use_real_fields so that we can switch back to real fields without initializing. And then this could be called by fields::change_m_number for m ≠ 1

else:
if self.m is not None:
self.fields.change_m_number(m)

def change_sources(self, new_sources):
"""
Change the list of sources in `Simulation.sources` to `new_sources`, and changes
Expand Down
9 changes: 9 additions & 0 deletions src/fields.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -793,4 +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) {
m = new_m;
for (int i = 0; i < num_chunks; i++) {
chunks[i]->change_m_number(new_m);
}
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

At the very least, this needs to check that new_m is consistent with whether we have real fields and abort if not. Better yet, have it call if (new_m != 0) use_real_fields(false); as noted above.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What do you envision use_real_fields(false) does in the case that the current simulation does have real fields (e.g. suppose originally m=0 but we are changing it to m=1). Wouldn't we need to reinitialize everything?

Or are you suggesting we refactor the initialization of the fields to use_real_fields(), such that

  • use_real_fields(true) deletes the extra array (the current behavior of use_real_fields())
  • use_real_fields(false) reallocates the fields if needed
  • does nothing if the simulation state is consistent

As a first pass, it might be easiest to leave use_real_fields as is and simply ensure there's a proper check in place for change_m_number() (like you suggest).

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My latest commit adds a check to make sure that the fields are consistent. If the user requests complex fields when the current setup is using real fields, I just abort.

This isn't a problem for the adjoint code. We have a similar check in the python routine. But rather than aborting, we just reinitialize.

}

void fields_chunk::change_m_number(double new_m) { m = new_m; }

} // namespace meep
2 changes: 2 additions & 0 deletions src/meep.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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<double> time_spent_on(time_sink sink);
Expand Down
35 changes: 17 additions & 18 deletions src/meepgeom.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2651,35 +2651,35 @@ 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);
LOOP_OVER_DIRECTIONS(dim, d) {
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;

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;
return fields_f * (row_1[dir_idx] - row_2[dir_idx]) / (2 * du);
}
// materials have some dispersion
else {
double orig = u[idx];
std::complex<double> row_1[3], row_2[3], dA_du[3];
Expand All @@ -2689,9 +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;

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);
}
}

Expand Down Expand Up @@ -2720,8 +2719,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) *
Expand Down