Check for invalid sources at r=0 in cylindrical coordinates (#2392)

* check for invalid sources at r=0 in cylindrical coordinates * exclude r=0 in source if |m| > 1 for test_adjoint_cyl.py * move check from add_volume_source to src_vol_chunkloop * Apply suggestions from code review * Update src/sources.cpp * fix formatting --------- Co-authored-by: Steven G. Johnson <[email protected]>
NanoComp · Feb 10, 2023 · 6abb8b3 · 6abb8b3
1 parent 395d821
commit 6abb8b3
Show file tree

Hide file tree

Showing 2 changed files with 43 additions and 6 deletions.
diff --git a/python/tests/test_adjoint_cyl.py b/python/tests/test_adjoint_cyl.py
@@ -36,10 +36,6 @@
 fcen = 1 / 1.55
 width = 0.2
 fwidth = width * fcen
-source_center = [design_r / 2, 0, -(sz / 2 - dpml + design_z / 2) / 2]
-source_size = mp.Vector3(design_r, 0, 0)
-src = mp.GaussianSource(frequency=fcen, fwidth=fwidth)
-source = [mp.Source(src, component=mp.Er, center=source_center, size=source_size)]
 
 ## random design region
 p = 0.5 * rng.rand(Nr * Nz)
@@ -48,6 +44,31 @@
 dp = deps * rng.rand(Nr * Nz)
 
 
+def get_source(m):
+    source_center = mp.Vector3(
+        design_r / 2,
+        0,
+        -(sz / 2 - dpml + design_z / 2) / 2,
+    )
+    # exclude r=0 whenever |m| > 1
+    source_size = mp.Vector3(
+        design_r if abs(m) <= 1 else 0.5 * design_r,
+        0,
+        0,
+    )
+    src = mp.GaussianSource(frequency=fcen, fwidth=fwidth)
+    source = [
+        mp.Source(
+            src,
+            component=mp.Er,
+            center=source_center,
+            size=source_size,
+        )
+    ]
+
+    return source
+
+
 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)
@@ -65,7 +86,7 @@ def forward_simulation(design_params, m, far_x):
         resolution=resolution,
         cell_size=cell_size,
         boundary_layers=boundary_layers,
-        sources=source,
+        sources=get_source(m),
         geometry=geometry,
         dimensions=dimensions,
         m=m,
@@ -111,7 +132,7 @@ def adjoint_solver(design_params, m, far_x):
         cell_size=cell_size,
         boundary_layers=boundary_layers,
         geometry=geometry,
-        sources=source,
+        sources=get_source(m),
         resolution=resolution,
         dimensions=dimensions,
         m=m,

diff --git a/src/sources.cpp b/src/sources.cpp
@@ -272,6 +272,22 @@ static void src_vol_chunkloop(fields_chunk *fc, int ichunk, component c, ivec is
     vec rel_loc = loc - data->center;
     amps_array[idx_vol] = IVEC_LOOP_WEIGHT(s0, s1, e0, e1, 1) * amp * data->A(rel_loc);
 
+    // check for invalid sources at r=0 in cylindrical coordinates
+    if (fc->gv.dim == Dcyl && loc.r() == 0 && amps_array[idx_vol] != 0.0) {
+      if (fc->m == 0 && (component_direction(c) == R || component_direction(c) == P))
+        meep::abort("Not possible to place a %s source at r=0 in "
+                    "cylindrical coordinates for m = 0.",
+                    component_name(c));
+      else if (fabs(fc->m) == 1.0 && component_direction(c) == Z)
+        meep::abort("Not possible to place a %s source at r=0 in "
+                    "cylindrical coordinates for |m| = 1.0.",
+                    component_name(c));
+      else
+        meep::abort("Not possible to place a source at r=0 in "
+                    "cylindrical coordinates for m = %g.",
+                    fc->m);
+    }
+
     /* for "D" sources, multiply by epsilon.  FIXME: this is not quite
        right because it doesn't handle non-diagonal chi1inv!
        similarly, for "B" sources, multiply by mu. */