From 6abb8b3c6692fb17d08a6a5b36b320e1c54f9c55 Mon Sep 17 00:00:00 2001
From: Ardavan Oskooi <ardavano@google.com>
Date: Thu, 9 Feb 2023 17:04:50 -0800
Subject: [PATCH] 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 <stevenj@mit.edu>
---
 python/tests/test_adjoint_cyl.py | 33 ++++++++++++++++++++++++++------
 src/sources.cpp                  | 16 ++++++++++++++++
 2 files changed, 43 insertions(+), 6 deletions(-)

diff --git a/python/tests/test_adjoint_cyl.py b/python/tests/test_adjoint_cyl.py
index 2ff6c0884..0ef340b41 100644
--- 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
index 4c872ba5f..50047bc31 100644
--- 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. */