Skip to content

Commit

Permalink
restrict weights parameter of MaterialGrid to range [0, 1] (#1881)
Browse files Browse the repository at this point in the history 
* restrict weights parameter of MaterialGrid to range [0,1]

* ensure that weights array of unit tests for adjoint solver have values less than 1

* add a check_weights method to MaterialGrid class
  • Loading branch information
@oskooi
oskooi authored Jan 5, 2022
1 parent 5dc6e68 commit 2e563e7
Show file tree
Hide file tree
Showing 5 changed files with 92 additions and 77 deletions.
61 changes: 31 additions & 30 deletions doc/docs/Python_User_Interface.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4395,10 +4395,9 @@ class MaterialGrid(object):

<div class="class_docstring" markdown="1">

This class is used to specify materials interpolated from discrete points on a rectilinear grid.
A class object is passed as the `material` argument of a [`Block`](#block) geometric object or
the `default_material` argument of the [`Simulation`](#Simulation) constructor (similar to a
[material function](#medium)).
This class is used to specify materials on a rectilinear grid. A class object is passed
as the `material` argument of a [`Block`](#block) geometric object or the `default_material`
argument of the [`Simulation`](#Simulation) constructor (similar to a [material function](#medium)).

</div>

Expand Down Expand Up @@ -4426,19 +4425,21 @@ def __init__(self,
Creates a `MaterialGrid` object.

The input are two materials `medium1` and `medium2` along with a weight function $u(x)$ which
is defined on discrete points of a rectilinear grid by the NumPy array `weights` of size `grid_size`
(a 3-tuple or `Vector3` of integers $N_x$,$N_y$,$N_z$). The resolution of the grid may be nonuniform
depending on the `size` property of the `Block` object as shown in the following example for a 2d
`MaterialGrid` with $N_x=5$ and $N_y=4$. $N_z=0$ implies that the `MaterialGrid` is extruded in the
$z$ direction. The grid points are defined at the corners of the voxels.
is defined on a rectilinear grid by the NumPy array `weights` of size `grid_size` (a 3-tuple or
`Vector3` of integers $N_x$,$N_y$,$N_z$). The resolution of the grid may be nonuniform depending
on the `size` property of the `Block` object as shown in the following example for a 2d `MaterialGrid`
with $N_x=5$ and $N_y=4$. $N_z=0$ implies that the `MaterialGrid` is extruded in the $z$ direction.
The grid points are defined at the corners of the voxels.

![](images/material_grid.png)

Elements of the `weights` array must be in the range [0,1] where 0 is `medium1` and 1 is `medium2`.
Two material types are supported: (1) frequency-independent isotropic $\varepsilon$ or $\mu$
and (2) `LorentzianSusceptibility`. `medium1` and `medium2` must both be the same type. The
materials are [bilinearly interpolated](https://en.wikipedia.org/wiki/Bilinear_interpolation)
from the rectilinear grid to Meep's [Yee grid](Yee_Lattice.md).
The `weights` array is used to define a linear interpolation from `medium1` to `medium2`.
Two material types are supported: (1) frequency-independent isotropic $\varepsilon$ (`epsilon_diag`
and `epsilon_offdiag` are interpolated) and (2) `LorentzianSusceptibility` (`sigma` and `sigma_offdiag`
are interpolated). `medium1` and `medium2` must both be the same type. The materials are
[bilinearly interpolated](https://en.wikipedia.org/wiki/Bilinear_interpolation) from the rectilinear
grid to Meep's [Yee grid](Yee_Lattice.md).

For improving accuracy, [subpixel smoothing](Subpixel_Smoothing.md) can be enabled by specifying
`do_averaging=True`. If you want to use a material grid to define a (nearly) discontinuous,
Expand All @@ -4455,9 +4456,9 @@ a *discontinuous* function from otherwise continuously varying (via the bilinear
grid values. Subpixel smoothing is fast and accurate because it exploits an analytic formulation
for level-set functions.

A nonzero damping term creates an artificial conductivity σ = u(1-u)*damping, which acts as
dissipation loss that penalize intermediate pixel values of non-binarized structures. The value of
damping should be proportional to 2π times the typical frequency of the problem.
A nonzero `damping` term creates an artificial conductivity $\sigma = u(1-u)*$`damping`, which acts as
dissipation loss that penalizes intermediate pixel values of non-binarized structures. The value of
`damping` should be proportional to $2\pi$ times the typical frequency of the problem.

It is possible to overlap any number of different `MaterialGrid`s. This can be useful for defining
grids which are symmetric (e.g., mirror, rotation). One way to set this up is by overlapping a
Expand Down Expand Up @@ -4499,7 +4500,7 @@ class Susceptibility(object):
<div class="class_docstring" markdown="1">

Parent class for various dispersive susceptibility terms, parameterized by an
anisotropic amplitude σ. See [Material Dispersion](Materials.md#material-dispersion).
anisotropic amplitude $\sigma$. See [Material Dispersion](Materials.md#material-dispersion).

</div>

Expand All @@ -4518,13 +4519,13 @@ def __init__(self,

<div class="method_docstring" markdown="1">

+ **`sigma` [`number`]** — The scale factor σ.
+ **`sigma` [`number`]** — The scale factor $\sigma$.

You can also specify an anisotropic σ tensor by using the property `sigma_diag`
which takes three numbers or a `Vector3` to give the σ$_n$ tensor diagonal, and
You can also specify an anisotropic $\sigma$ tensor by using the property `sigma_diag`
which takes three numbers or a `Vector3` to give the $\sigma_n$ tensor diagonal, and
`sigma_offdiag` which specifies the offdiagonal elements (defaults to 0). That is,
`sigma_diag=mp.Vector3(a, b, c)` and `sigma_offdiag=mp.Vector3(u, v, w)`
corresponds to a σ tensor
corresponds to a $\sigma$ tensor

\begin{pmatrix} a & u & v \\ u & b & w \\ v & w & c \end{pmatrix}

Expand All @@ -4545,7 +4546,7 @@ class LorentzianSusceptibility(Susceptibility):

Specifies a single dispersive susceptibility of Lorentzian (damped harmonic
oscillator) form. See [Material Dispersion](Materials.md#material-dispersion), with
the parameters (in addition to σ):
the parameters (in addition to $\sigma$):

</div>

Expand All @@ -4563,7 +4564,7 @@ def __init__(self, frequency=0.0, gamma=0.0, **kwargs):

+ **`frequency` [`number`]** — The resonance frequency $f_n = \omega_n / 2\pi$.

+ **`gamma` [`number`]** — The resonance loss rate $γ_n / 2\pi$.
+ **`gamma` [`number`]** — The resonance loss rate $\gamma_n / 2\pi$.

Note: multiple objects with identical values for the `frequency` and `gamma` but
different `sigma` will appear as a *single* Lorentzian susceptibility term in the
Expand All @@ -4585,7 +4586,7 @@ class DrudeSusceptibility(Susceptibility):
<div class="class_docstring" markdown="1">

Specifies a single dispersive susceptibility of Drude form. See [Material
Dispersion](Materials.md#material-dispersion), with the parameters (in addition to σ):
Dispersion](Materials.md#material-dispersion), with the parameters (in addition to $\sigma$):

</div>

Expand All @@ -4602,9 +4603,9 @@ def __init__(self, frequency=0.0, gamma=0.0, **kwargs):
<div class="method_docstring" markdown="1">

+ **`frequency` [`number`]** — The frequency scale factor $f_n = \omega_n / 2\pi$
which multiplies σ (not a resonance frequency).
which multiplies $\sigma$ (not a resonance frequency).

+ **`gamma` [`number`]** — The loss rate $γ_n / 2\pi$.
+ **`gamma` [`number`]** — The loss rate $\gamma_n / 2\pi$.

</div>

Expand Down Expand Up @@ -4645,13 +4646,13 @@ def __init__(self,

<div class="method_docstring" markdown="1">

+ **`sigma` [`number`]** — The scale factor σ.
+ **`sigma` [`number`]** — The scale factor $\sigma$.

You can also specify an anisotropic σ tensor by using the property `sigma_diag`
which takes three numbers or a `Vector3` to give the σ$_n$ tensor diagonal, and
You can also specify an anisotropic $\sigma$ tensor by using the property `sigma_diag`
which takes three numbers or a `Vector3` to give the $\sigma_n$ tensor diagonal, and
`sigma_offdiag` which specifies the offdiagonal elements (defaults to 0). That is,
`sigma_diag=mp.Vector3(a, b, c)` and `sigma_offdiag=mp.Vector3(u, v, w)`
corresponds to a σ tensor
corresponds to a $\sigma$ tensor

\begin{pmatrix} a & u & v \\ u & b & w \\ v & w & c \end{pmatrix}

Expand Down
78 changes: 48 additions & 30 deletions python/geom.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -524,29 +524,46 @@ def _get_epsmu(self, diag, offdiag, susceptibilities, conductivity_diag, conduct

class MaterialGrid(object):
"""
This class is used to specify materials interpolated from discrete points on a rectilinear grid.
A class object is passed as the `material` argument of a [`Block`](#block) geometric object or
the `default_material` argument of the [`Simulation`](#Simulation) constructor (similar to a
[material function](#medium)).
This class is used to specify materials on a rectilinear grid. A class object is passed
as the `material` argument of a [`Block`](#block) geometric object or the `default_material`
argument of the [`Simulation`](#Simulation) constructor (similar to a [material function](#medium)).
"""
def __init__(self,grid_size,medium1,medium2,weights=None,grid_type="U_DEFAULT",do_averaging=True,beta=0,eta=0.5,damping=0):
def check_weights(self, w):
if np.amin(w) < 0. or np.amax(w) > 1.:
warnings.warn('The weights parameter of MaterialGrid must be in the range [0,1].')
return np.clip(w, 0., 1.)
else:
return w

def __init__(self,
grid_size,
medium1,
medium2,
weights=None,
grid_type="U_DEFAULT",
do_averaging=True,
beta=0,
eta=0.5,
damping=0):
"""
Creates a `MaterialGrid` object.
The input are two materials `medium1` and `medium2` along with a weight function $u(x)$ which
is defined on discrete points of a rectilinear grid by the NumPy array `weights` of size `grid_size`
(a 3-tuple or `Vector3` of integers $N_x$,$N_y$,$N_z$). The resolution of the grid may be nonuniform
depending on the `size` property of the `Block` object as shown in the following example for a 2d
`MaterialGrid` with $N_x=5$ and $N_y=4$. $N_z=0$ implies that the `MaterialGrid` is extruded in the
$z$ direction. The grid points are defined at the corners of the voxels.
is defined on a rectilinear grid by the NumPy array `weights` of size `grid_size` (a 3-tuple or
`Vector3` of integers $N_x$,$N_y$,$N_z$). The resolution of the grid may be nonuniform depending
on the `size` property of the `Block` object as shown in the following example for a 2d `MaterialGrid`
with $N_x=5$ and $N_y=4$. $N_z=0$ implies that the `MaterialGrid` is extruded in the $z$ direction.
The grid points are defined at the corners of the voxels.
![](images/material_grid.png)
Elements of the `weights` array must be in the range [0,1] where 0 is `medium1` and 1 is `medium2`.
Two material types are supported: (1) frequency-independent isotropic $\\varepsilon$ or $\\mu$
and (2) `LorentzianSusceptibility`. `medium1` and `medium2` must both be the same type. The
materials are [bilinearly interpolated](https://en.wikipedia.org/wiki/Bilinear_interpolation)
from the rectilinear grid to Meep's [Yee grid](Yee_Lattice.md).
The `weights` array is used to define a linear interpolation from `medium1` to `medium2`.
Two material types are supported: (1) frequency-independent isotropic $\\varepsilon$ (`epsilon_diag`
and `epsilon_offdiag` are interpolated) and (2) `LorentzianSusceptibility` (`sigma` and `sigma_offdiag`
are interpolated). `medium1` and `medium2` must both be the same type. The materials are
[bilinearly interpolated](https://en.wikipedia.org/wiki/Bilinear_interpolation) from the rectilinear
grid to Meep's [Yee grid](Yee_Lattice.md).
For improving accuracy, [subpixel smoothing](Subpixel_Smoothing.md) can be enabled by specifying
`do_averaging=True`. If you want to use a material grid to define a (nearly) discontinuous,
Expand All @@ -563,9 +580,9 @@ def __init__(self,grid_size,medium1,medium2,weights=None,grid_type="U_DEFAULT",d
grid values. Subpixel smoothing is fast and accurate because it exploits an analytic formulation
for level-set functions.
A nonzero damping term creates an artificial conductivity σ = u(1-u)*damping, which acts as
dissipation loss that penalize intermediate pixel values of non-binarized structures. The value of
damping should be proportional to 2π times the typical frequency of the problem.
A nonzero `damping` term creates an artificial conductivity $\\sigma = u(1-u)*$`damping`, which acts as
dissipation loss that penalizes intermediate pixel values of non-binarized structures. The value of
`damping` should be proportional to $2\\pi$ times the typical frequency of the problem.
It is possible to overlap any number of different `MaterialGrid`s. This can be useful for defining
grids which are symmetric (e.g., mirror, rotation). One way to set this up is by overlapping a
Expand Down Expand Up @@ -596,7 +613,7 @@ def isclose(a, b, rel_tol=1e-09, abs_tol=0.0):
elif weights.size != self.num_params:
raise ValueError("weights of shape {} do not match user specified grid dimension: {}".format(weights.size,self.grid_size))
else:
self.weights = weights.flatten().astype(np.float64)
self.weights = self.check_weights(weights).flatten().astype(np.float64)

grid_type_dict = {
"U_MIN":0,
Expand All @@ -609,28 +626,29 @@ def isclose(a, b, rel_tol=1e-09, abs_tol=0.0):
self.grid_type = grid_type_dict[grid_type]

self.swigobj = None
def update_weights(self,x):

def update_weights(self, x):
"""
Reset the `weights` to `x`.
"""
if x.size != self.num_params:
raise ValueError("weights of shape {} do not match user specified grid dimension: {}".format(self.weights.size,self.grid_size))
self.weights[:]=x.flatten().astype(np.float64)
self.weights[:] = self.check_weights(x).flatten().astype(np.float64)

class Susceptibility(object):
"""
Parent class for various dispersive susceptibility terms, parameterized by an
anisotropic amplitude σ. See [Material Dispersion](Materials.md#material-dispersion).
anisotropic amplitude $\\sigma$. See [Material Dispersion](Materials.md#material-dispersion).
"""
def __init__(self, sigma_diag=Vector3(), sigma_offdiag=Vector3(), sigma=None):
"""
+ **`sigma` [`number`]** — The scale factor σ.
+ **`sigma` [`number`]** — The scale factor $\\sigma$.
You can also specify an anisotropic σ tensor by using the property `sigma_diag`
which takes three numbers or a `Vector3` to give the σ$_n$ tensor diagonal, and
You can also specify an anisotropic $\\sigma$ tensor by using the property `sigma_diag`
which takes three numbers or a `Vector3` to give the $\\sigma_n$ tensor diagonal, and
`sigma_offdiag` which specifies the offdiagonal elements (defaults to 0). That is,
`sigma_diag=mp.Vector3(a, b, c)` and `sigma_offdiag=mp.Vector3(u, v, w)`
corresponds to a σ tensor
corresponds to a $\\sigma$ tensor
\\begin{pmatrix} a & u & v \\\\ u & b & w \\\\ v & w & c \\end{pmatrix}
"""
Expand All @@ -648,13 +666,13 @@ class LorentzianSusceptibility(Susceptibility):
"""
Specifies a single dispersive susceptibility of Lorentzian (damped harmonic
oscillator) form. See [Material Dispersion](Materials.md#material-dispersion), with
the parameters (in addition to σ):
the parameters (in addition to $\\sigma$):
"""
def __init__(self, frequency=0.0, gamma=0.0, **kwargs):
"""
+ **`frequency` [`number`]** — The resonance frequency $f_n = \\omega_n / 2\\pi$.
+ **`gamma` [`number`]** — The resonance loss rate $γ_n / 2\\pi$.
+ **`gamma` [`number`]** — The resonance loss rate $\\gamma_n / 2\\pi$.
Note: multiple objects with identical values for the `frequency` and `gamma` but
different `sigma` will appear as a *single* Lorentzian susceptibility term in the
Expand All @@ -675,14 +693,14 @@ def eval_susceptibility(self,freq):
class DrudeSusceptibility(Susceptibility):
"""
Specifies a single dispersive susceptibility of Drude form. See [Material
Dispersion](Materials.md#material-dispersion), with the parameters (in addition to σ):
Dispersion](Materials.md#material-dispersion), with the parameters (in addition to $\\sigma$):
"""
def __init__(self, frequency=0.0, gamma=0.0, **kwargs):
"""
+ **`frequency` [`number`]** — The frequency scale factor $f_n = \\omega_n / 2\\pi$
which multiplies σ (not a resonance frequency).
which multiplies $\\sigma$ (not a resonance frequency).
+ **`gamma` [`number`]** — The loss rate $γ_n / 2\\pi$.
+ **`gamma` [`number`]** — The loss rate $\\gamma_n / 2\\pi$.
"""
super(DrudeSusceptibility, self).__init__(**kwargs)
self.frequency = frequency
Expand Down
2 changes: 1 addition & 1 deletion python/tests/test_adjoint_cyl.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@
size=source_size)]

## random design region
p = rng.rand(Nr*Nz)
p = 0.5*rng.rand(Nr*Nz)
## random epsilon perturbation for design region
deps = 1e-5
dp = deps*rng.rand(Nr*Nz)
Expand Down
12 changes: 6 additions & 6 deletions python/tests/test_adjoint_jax.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -182,17 +182,17 @@ def test_design_region_monitor_helpers(self):
class WrapperTest(ApproxComparisonTestCase):
@parameterized.parameterized.expand([
('1500_1550bw_01relative_gaussian_port1',
onp.linspace(1 / 1.50, 1 / 1.55, 3).tolist(), 0.1, 1.0, 0),
onp.linspace(1 / 1.50, 1 / 1.55, 3).tolist(), 0.1, 0.5, 0),
('1550_1600bw_02relative_gaussian_port1',
onp.linspace(1 / 1.55, 1 / 1.60, 3).tolist(), 0.2, 1.0, 0),
onp.linspace(1 / 1.55, 1 / 1.60, 3).tolist(), 0.2, 0.5, 0),
('1500_1600bw_03relative_gaussian_port1',
onp.linspace(1 / 1.50, 1 / 1.60, 4).tolist(), 0.3, 1.0, 0),
onp.linspace(1 / 1.50, 1 / 1.60, 4).tolist(), 0.3, 0.5, 0),
('1500_1550bw_01relative_gaussian_port2',
onp.linspace(1 / 1.50, 1 / 1.55, 3).tolist(), 0.1, 1.0, 1),
onp.linspace(1 / 1.50, 1 / 1.55, 3).tolist(), 0.1, 0.5, 1),
('1550_1600bw_02relative_gaussian_port2',
onp.linspace(1 / 1.55, 1 / 1.60, 3).tolist(), 0.2, 1.0, 1),
onp.linspace(1 / 1.55, 1 / 1.60, 3).tolist(), 0.2, 0.5, 1),
('1500_1600bw_03relative_gaussian_port2',
onp.linspace(1 / 1.50, 1 / 1.60, 4).tolist(), 0.3, 1.0, 1),
onp.linspace(1 / 1.50, 1 / 1.60, 4).tolist(), 0.3, 0.5, 1),
])
def test_wrapper_gradients(self, _, frequencies, gaussian_rel_width,
design_variable_fill_value, excite_port_idx):
Expand Down
Loading

0 comments on commit 2e563e7

Please sign in to comment.