search

NanoComp / meep

free finite-difference time-domain (FDTD) software for electromagnetic simulations
GNU General Public License v2.0
1.25k stars 626 forks source link

fields instability when combining 2+ mirror symmetries with periodic boundary conditions in 3d #132

Open oskooi opened 7 years ago

oskooi commented 7 years ago

As has been already reported on the meep-discuss list, there is a bug which leads to field instabilities when combining 3 mirror-symmetry objects with periodic boundary conditions in 3d.

This is demonstrated by the following simple example of a point source in vacuum:

(set-param! resolution 10)
(set! geometry-lattice (make lattice (size 1 1 1)))

(set! sources (list (src (make source (src (make gaussian-src (frequency 1) (fwidth 1))) (component Ez) (center 0 0 0)))))

(set! k-point (vector3 0 0 0))

(set! symmetries (list (make mirror-sym (direction X))
                       (make mirror-sym (direction Y))
                       (make mirror-sym (direction Z) (phase -1))))

(define print-field (lambda () (print "ez:, " (meep-time) ", " (magnitude (get-field-point Ez (vector3 0 0 0))) "\n")))

(run-until 500 (at-every 3.4 print-field))

For this script, the fields begin to blow up early on:

-----------
Initializing structure...
Working in 3D dimensions.
Computational cell is 1 x 1 x 1 with resolution 10
Halving computational cell along direction x
Halving computational cell along direction y
Halving computational cell along direction z
time for set_epsilon = 0.000983953 s
-----------
ez:, 3.4000000000000004, 0.11299371272955527
ez:, 6.800000000000001, 1.3292101237116931
ez:, 10.200000000000001, 2.402271717485686
ez:, 13.600000000000001, 4.484463313924006
ez:, 17.0, 3.8996648343918103
ez:, 20.400000000000002, 0.13634058216335865
ez:, 23.8, 4.629265409516416
ez:, 27.200000000000003, 4.998705102736996
ez:, 30.6, 1.2298673000561804
ez:, 34.0, 2.6981516161770243
ez:, 37.4, 4.377899104931889
ez:, 40.800000000000004, 3.5929413491696733
ez:, 44.2, 15.319503849169674
ez:, 47.6, 1194.3195038491704
ez:, 51.0, 2915498.319503849
ez:, 54.400000000000006, 559645525.6804962
ez:, 57.800000000000004, 4232176173909.6807

There are no instabilities if the Z mirror-symmetry object is removed:

-----------
Initializing structure...
Working in 3D dimensions.
Computational cell is 1 x 1 x 1 with resolution 10
Halving computational cell along direction x
Halving computational cell along direction y
time for set_epsilon = 0.00132418 s
-----------
ez:, 3.4000000000000004, 0.11299371272955483
ez:, 6.800000000000001, 1.3292101237116927
ez:, 10.200000000000001, 2.402271717485691
ez:, 13.600000000000001, 4.484463313924063
ez:, 17.0, 3.899664834392148
ez:, 20.400000000000002, 0.13634058216075873
ez:, 23.8, 4.629265409492891
ez:, 27.200000000000003, 4.998705102472007
ez:, 30.6, 1.2298672971468396
ez:, 34.0, 2.6981515888410907
ez:, 37.4, 4.3779221161355615
ez:, 40.800000000000004, 3.6004200025033697
ez:, 44.2, 0.1855690755252214
ez:, 47.6, 4.1906787491702815
ez:, 51.0, 5.302888125953832
ez:, 54.400000000000006, 1.6628091344931328
ez:, 57.800000000000004, 2.9025704016780374
ChristopherHogan commented 6 years ago

Triggering the bug requires periodic boundary conditions, and more than one mirror symmetry, one of which must be Z.

Symmetry Output
[Mirror(X), Mirror(Y)] correct
[Mirror(Z)] correct
[Mirror(X), Mirror(Z)] incorrect
[Mirror(Y), Mirror(Z)] incorrect
[Mirror(X), Mirror(Y), Mirror(Z)] incorrect

Also, based on the formulas Steven gave me for converting from 3D points to a 1D array, there appears to be an off-by-one error. For example, with only X and Y mirror symmetries, fields::chunks[0]->sources[2].index equals {598, 599}. Converting those indices to 3D points:

# gv size (8, 8, 11)
# 598
iz = 598 % 11 = 4
iy = (598 / 11) % 8 = 54 % 8 = 6
ix = (598 / 11) / 8 = 54 / 8 = 6
=> (6, 6, 4)

# 599
iz = 599 % 11 = 5
iy = (599 / 11) % 8 = 54 % 8 = 6
ix = (599 / 11) / 8 = 54 / 8 = 6
=> (6, 6, 5)

Now converting the 3D points into a 1D array with Z symmetry:

# gv size (8, 8, 8)
(6*8 + 6)*8 + 4 = 54*8 + 4 = 436
(6*8 + 6)*8 + 5 = 54*8 + 5 = 437

However, adding Z symmetry, fields::chunks[0]->sources[2].index equals {437, 438}, not {436, 437}. 437 and 438 correspond to (6, 6, 5) and (6, 6, 6).

stevengj commented 5 years ago

By any chance, does this patch fix it?

--- a/src/structure.cpp
+++ b/src/structure.cpp
@@ -160,7 +160,7 @@ void structure::choose_chunkdivision(const grid_volume &thegv, int desired_num_c
     v = gv.surroundings();
     // Pad the little cell in any direction that we've shrunk:
     for (int d = 0; d < 3; d++)
-      if (break_this[d]) gv = gv.pad((direction)d);
+      if (break_this[d] && thegv.num_direction((direction)d) & 1) gv = gv.pad((direction)d);
   }
stevengj commented 5 years ago

(In #1045, @oskooi said that the problem arises when the Padding %s to even number of grid points message is not printed, so it makes me wonder whether the problem has to do with the gv.pad call above, or maybe something in gv.halve for even numbers of grid points.)

Another thing to check would be this line: it says note that icenter-io is always even by construction of grid_volume::icenter. Would be good to check that this is actually true by adding an assertion or something on that line.

(That comment goes way back to commit 4e8a6d51aa104f7bfefb74242923455b16137a3d … in @droundy's original implementation, symmetries were only supported for an even number of pixels in a given direction. I "fixed" this, but maybe I messed something up at that time?)

oskooi commented 5 years ago

Adding the line master_printf("icenter-io: %d\n",icenter().in_direction(d) - io.in_direction(d)); (without the patch above) just before src/vec.cpp:1150 seems to always show even output for the test described in #1043:comment.

resolution = 15

Using MPI version 3.1, 3 processes
-----------
Initializing structure...
Padding y to even number of grid points.
Padding z to even number of grid points.
Halving computational cell along direction y
icenter-io: 76
Halving computational cell along direction z
icenter-io: 76
Splitting into 3 chunks evenly
time for choose_chunkdivision = 0.000445383 s
Working in 3D dimensions.
Computational cell is 7 x 5 x 5 with resolution 15
time for set_epsilon = 0.0418544 s
-----------
Meep progress: 62.53333333333333/70.0 = 89.3% done in 4.0s, 0.5s to go
on time step 1876 (time=62.5333), 0.00213281 s/step
run 0 finished at t = 70.0 (2100 timesteps)
flux:, 1.522463, 1.522463, -0.000000, 0.000000, -0.000000, 0.000000

resolution = 20 (fields blow up)

sing MPI version 3.1, 3 processes
-----------
Initializing structure...
Halving computational cell along direction y
icenter-io: 100
Halving computational cell along direction z
icenter-io: 100
Splitting into 3 chunks evenly
time for choose_chunkdivision = 0.00058396 s
Working in 3D dimensions.
Computational cell is 7 x 5 x 5 with resolution 20
time for set_epsilon = 0.0899128 s
-----------
Meep progress: 18.05/70.0 = 25.8% done in 4.0s, 11.5s to go
on time step 722 (time=18.05), 0.0055451 s/step
Meep progress: 38.225/70.0 = 54.6% done in 8.0s, 6.7s to go
on time step 1529 (time=38.225), 0.00496038 s/step
Meep progress: 58.425000000000004/70.0 = 83.5% done in 12.0s, 2.4s to go
on time step 2337 (time=58.425), 0.00495388 s/step
run 0 finished at t = 70.0 (2800 timesteps)
flux:, 1820348954164070400.000000, 1081458856371996928.000000, -82773801567460112.000000, 82773801567460112.000000, 8971776376992865280.000000, -8971776376992865280.000000

resolution = 30 (fields blow up)

Using MPI version 3.1, 3 processes
-----------
Initializing structure...
Halving computational cell along direction y
icenter-io: 150
Halving computational cell along direction z
icenter-io: 150
Splitting into 3 chunks evenly
time for choose_chunkdivision = 0.000662629 s
Working in 3D dimensions.
Computational cell is 7 x 5 x 5 with resolution 30
time for set_epsilon = 0.149383 s
-----------
Meep progress: 3.0833333333333335/70.0 = 4.4% done in 4.0s, 87.1s to go
on time step 185 (time=3.08333), 0.0217037 s/step
Meep progress: 7.1499999999999995/70.0 = 10.2% done in 8.0s, 70.6s to go
on time step 429 (time=7.15), 0.0164605 s/step
Meep progress: 11.2/70.0 = 16.0% done in 12.0s, 63.2s to go
on time step 672 (time=11.2), 0.0165022 s/step
Meep progress: 15.299999999999999/70.0 = 21.9% done in 16.0s, 57.4s to go
on time step 918 (time=15.3), 0.01629 s/step
Meep progress: 19.316666666666666/70.0 = 27.6% done in 20.1s, 52.6s to go
on time step 1159 (time=19.3167), 0.0166254 s/step
Meep progress: 23.266666666666666/70.0 = 33.2% done in 24.1s, 48.3s to go
on time step 1396 (time=23.2667), 0.0169003 s/step
Meep progress: 27.28333333333333/70.0 = 39.0% done in 28.1s, 43.9s to go
on time step 1637 (time=27.2833), 0.016599 s/step
Meep progress: 31.25/70.0 = 44.6% done in 32.1s, 39.8s to go
on time step 1875 (time=31.25), 0.0168311 s/step
Meep progress: 35.266666666666666/70.0 = 50.4% done in 36.1s, 35.5s to go
on time step 2116 (time=35.2667), 0.0166289 s/step
Meep progress: 39.3/70.0 = 56.1% done in 40.1s, 31.3s to go
on time step 2358 (time=39.3), 0.0165408 s/step
Meep progress: 43.38333333333333/70.0 = 62.0% done in 44.1s, 27.0s to go
on time step 2603 (time=43.3833), 0.0163723 s/step
Meep progress: 47.46666666666667/70.0 = 67.8% done in 48.1s, 22.8s to go
on time step 2848 (time=47.4667), 0.0163545 s/step
Meep progress: 51.43333333333333/70.0 = 73.5% done in 52.1s, 18.8s to go
on time step 3086 (time=51.4333), 0.0168205 s/step
Meep progress: 55.46666666666667/70.0 = 79.2% done in 56.1s, 14.7s to go
on time step 3328 (time=55.4667), 0.0165867 s/step
Meep progress: 59.53333333333333/70.0 = 85.0% done in 60.1s, 10.6s to go
on time step 3572 (time=59.5333), 0.0164323 s/step
Meep progress: 63.6/70.0 = 90.9% done in 64.1s, 6.5s to go
on time step 3816 (time=63.6), 0.0164278 s/step
Meep progress: 67.63333333333333/70.0 = 96.6% done in 68.1s, 2.4s to go
on time step 4058 (time=67.6333), 0.0165797 s/step
run 0 finished at t = 70.0 (4200 timesteps)
flux:, 1459390448574672659383315722416858030678581706752.000000, 88649338771690060606454862583248934889125314560.000000, 174732717817583568259904759287750632314966638592.000000, -174732717817583568259904759287750632314966638592.000000, 1524448985721283365421259838898912231715185360896.000000, -1524448985721283365421259838898912231715185360896.000000
oskooi commented 4 years ago

I think I have discovered the cause of this bug as well as a potential fix/workaround.

First is the cause. The fields blow up whenever three conditions are satisfied: (1) the k_point is not False (i.e., Bloch-periodic boundary conditions), (2) there are at least two mirror-symmetry planes, and (3) one of the symmetries has odd phase. The third criteria is important and points the way to a fix.

(The original example in this thread involved an Ez source with an odd mirror-symmetry plane in z. The same bug can also be demonstrated using an Ex or Ey source with an odd mirror-symmetry in x or y, respectively.)

The fix: force cell padding in only the direction of the odd mirror-symmetry plane by increasing/decreasing the size of the cell by one pixel (i.e., for an Ej source in a cell with even number of pixels in all directions, add one pixel to the j direction). The key is to modify the cell size so that padding occurs in only the direction of the odd mirror-symmetry plane.

This is demonstrated by the following example for an Ey source with the fix shown in Case 2.

(set-param! resolution 10)

(define-param sx 1.0)
(define-param sy 1.0)
(define-param sz 1.0)
(set! geometry-lattice (make lattice (size sx sy sz)))

(set! sources (list (src (make source
                           (src (make gaussian-src (frequency 1) (fwidth 1)))
                           (component Ey)
                           (center 0 0 0)))))

(set! k-point (vector3 0 0 0))

(set! symmetries (list (make mirror-sym (direction X) (phase +1))
                       (make mirror-sym (direction Y) (phase -1))
                       (make mirror-sym (direction Z) (phase +1))))

(define print-field (lambda ()
                      (print "ey:, "
                             (meep-time) ", "
                             (magnitude (get-field-point Ey (vector3 0 0 0))) "\n")))

(run-until 100 (at-every 3.4 print-field))

Case 1: no padding

> meep sx=1.0 sy=1.0 sz=1.0 unstable.ctl
Using MPI version 3.1, 1 processes
command-line param: sx=1.0
command-line param: sy=1.0
command-line param: sz=1.0
-----------
Initializing structure...
Working in 3D dimensions.
Computational cell is 1 x 1 x 1 with resolution 10
Halving computational cell along direction x
Halving computational cell along direction y
Halving computational cell along direction z
time for choose_chunkdivision = 2.278e-05 s
time for set_epsilon = 0.000681023 s
-----------
..........
ey:, 91.80000000000001, 5.731901655445799e40
ey:, 95.2, 4.529551484666823e43
ey:, 98.60000000000001, 5.542780816649825e46
run 0 finished at t = 100.0 (2000 timesteps)

result: unstable

*** Case 2: padding in y

> meep sx=1.0 sy=1.1 sz=1.0 unstable.ctl
Using MPI version 3.1, 1 processes
command-line param: sx=1.0
command-line param: sy=1.1
command-line param: sz=1.0
-----------
Initializing structure...
Working in 3D dimensions.
Computational cell is 1 x 1.1 x 1 with resolution 10
Padding y to even number of grid points.
Halving computational cell along direction x
Halving computational cell along direction y
Halving computational cell along direction z
time for choose_chunkdivision = 1.9466e-05 s
time for set_epsilon = 0.000759672 s
-----------
..........
ey:, 91.80000000000001, 0.62514828007822
ey:, 95.2, 3.9391504102184896
ey:, 98.60000000000001, 5.080533636288392
run 0 finished at t = 100.0 (2000 timesteps)

*** result: stable

Case 3: padding in z

> meep sx=1.0 sy=1.0 sz=1.1 unstable.ctl
Using MPI version 3.1, 1 processes
command-line param: sx=1.0
command-line param: sy=1.0
command-line param: sz=1.1
-----------
Initializing structure...
Working in 3D dimensions.
Computational cell is 1 x 1 x 1.1 with resolution 10
Padding z to even number of grid points.
Halving computational cell along direction x
Halving computational cell along direction y
Halving computational cell along direction z
time for choose_chunkdivision = 2.2074e-05 s
time for set_epsilon = 0.000752958 s
-----------
..........
ey:, 91.80000000000001, 4.275517439564672e36
ey:, 95.2, 1.903504180320423e40
ey:, 98.60000000000001, 3.213799530049778e42
run 0 finished at t = 100.0 (2000 timesteps)

result: unstable

Case 4: padding in x

> meep sx=1.1 sy=1.0 sz=1.0 unstable.ctl
Using MPI version 3.1, 1 processes
command-line param: sx=1.1
command-line param: sy=1.0
command-line param: sz=1.0
-----------
Initializing structure...
Working in 3D dimensions.
Computational cell is 1.1 x 1 x 1 with resolution 10
Padding x to even number of grid points.
Halving computational cell along direction x
Halving computational cell along direction y
Halving computational cell along direction z
time for choose_chunkdivision = 2.1227e-05 s
time for set_epsilon = 0.000754938 s
-----------
..........
ey:, 91.80000000000001, 4.2233036407185267e37
ey:, 95.2, 8.552081671358194e39
ey:, 98.60000000000001, 6.986159356210885e41
run 0 finished at t = 100.0 (2000 timesteps)

result: unstable

Case 5: padding in x and y

> meep sx=1.1 sy=1.1 sz=1.0 unstable.ctl
Using MPI version 3.1, 1 processes
command-line param: sx=1.1
command-line param: sy=1.1
command-line param: sz=1.0
-----------
Initializing structure...
Working in 3D dimensions.
Computational cell is 1.1 x 1.1 x 1 with resolution 10
Padding x to even number of grid points.
Padding y to even number of grid points.
Halving computational cell along direction x
Halving computational cell along direction y
Halving computational cell along direction z
time for choose_chunkdivision = 2.2002e-05 s
time for set_epsilon = 0.000850702 s
-----------
..........
ey:, 91.80000000000001, 1.6522542159769548e25
ey:, 95.2, 7.130121909132e26
ey:, 98.60000000000001, 3.0769259323102425e28
run 0 finished at t = 100.0 (2000 timesteps)

result: unstable

Case 6: padding in y and z

> meep sx=1.0 sy=1.1 sz=1.1 unstable.ctl
Using MPI version 3.1, 1 processes
command-line param: sx=1.0
command-line param: sy=1.1
command-line param: sz=1.1
-----------
Initializing structure...
Working in 3D dimensions.
Computational cell is 1 x 1.1 x 1.1 with resolution 10
Padding y to even number of grid points.
Padding z to even number of grid points.
Halving computational cell along direction x
Halving computational cell along direction y
Halving computational cell along direction z
time for choose_chunkdivision = 2.5274e-05 s
time for set_epsilon = 0.000847353 s
-----------
..........
ey:, 91.80000000000001, 9.383694164491804e24
ey:, 95.2, 4.0494303300219355e26
ey:, 98.60000000000001, 1.747487259308995e28
run 0 finished at t = 100.0 (2000 timesteps)

result: unstable

Case 6: padding in x, y, and z

> meep sx=1.1 sy=1.1 sz=1.1 unstable.ctl
Using MPI version 3.1, 1 processes
command-line param: sx=1.1
command-line param: sy=1.1
command-line param: sz=1.1
-----------
Initializing structure...
Working in 3D dimensions.
Computational cell is 1.1 x 1.1 x 1.1 with resolution 10
Padding x to even number of grid points.
Padding y to even number of grid points.
Padding z to even number of grid points.
Halving computational cell along direction x
Halving computational cell along direction y
Halving computational cell along direction z
time for choose_chunkdivision = 2.5263e-05 s
time for set_epsilon = 0.000953828 s
-----------
..........
ey:, 91.80000000000001, 34272704218962.453
ey:, 95.2, 5078456364939438.0
ey:, 98.60000000000001, 191907472156637380.0
run 0 finished at t = 100.0 (2000 timesteps)

result: unstable

oskooi commented 4 years ago

With k_point of (0.5,0.5,0.5) the fields are unstable regardless of whether any or all directions padded.

(set-param! resolution 10)

(define-param sx 1.0)
(define-param sy 1.0)
(define-param sz 1.0)
(set! geometry-lattice (make lattice (size sx sy sz)))

(set! sources (list (src (make source
                           (src (make gaussian-src (frequency 1) (fwidth 1)))
                           (component Ey)
                           (center 0 0 0)))))

(set! k-point (vector3 0.5 0.5 0.5))                                                 

(set! symmetries (list (make mirror-sym (direction X) (phase +1))
                       (make mirror-sym (direction Y) (phase -1))
                       (make mirror-sym (direction Z) (phase +1))))

(define print-field (lambda ()
                      (print "ey:, "
                             (meep-time) ", "
                             (magnitude (get-field-point Ey (vector3 0 0 0))) "\n")))

(run-until 100 (at-every 3.4 print-field))
stevengj commented 4 years ago

Might be worth looking at the fields at the edge of the cell (using get_field_point and plugging in a Yee grid point).

(Note that if you call the integrate_fields function with a single component, it will call your function at the Yee grid points for that component.)