-
Notifications
You must be signed in to change notification settings - Fork 1
/
mesti_main.jl
1882 lines (1771 loc) · 115 KB
/
mesti_main.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# Export composite data types
export Source_struct
export Syst
# Export a function mesti()
export mesti
mutable struct Source_struct
# A composite data type to specfiy source
# See also: mesti and mesti2s
pos::Vector{Vector{Int64}}
data::Union{Vector{Array{Int64,2}},Vector{Array{Float64,2}},Vector{Array{ComplexF64,2}},
Vector{Array{Int64,3}},Vector{Array{Float64,3}},Vector{Array{ComplexF64,3}},
Vector{Array{Int64,4}},Vector{Array{Float64,4}},Vector{Array{ComplexF64,4}}}
ind::Vector{Vector{Int64}}
isempty::Integer
Source_struct() = new()
end
mutable struct Syst
# A composite data type to specfiy system
# See also: mesti and mesti2s
# Below are used in both mesti() and mesti2s()
epsilon_xx::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Matrix{Int64},Matrix{Float64},Matrix{ComplexF64},Nothing}
epsilon_xy::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}
epsilon_xz::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}
epsilon_yx::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}
epsilon_yy::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}
epsilon_yz::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}
epsilon_zx::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}
epsilon_zy::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}
epsilon_zz::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}
length_unit::String
wavelength::Number
dx::Real
xBC::Union{String,Nothing}
yBC::String
kx_B::Union{Number,Nothing}
ky_B::Number
# Below are used in mesti() only
zBC::String
kz_B::Number
PML::Union{PML,Vector{PML}}
PML_type::String
# Below are used in mesti2s() only
epsilon_low::Union{Real,Nothing}
epsilon_high::Union{Real,Nothing}
zPML::Union{PML,Vector{PML},Nothing}
Syst() = new()
end
"""
MESTI Multi-source frequency-domain electromagnetic simulations.
---3D field profile---
(Ex, Ey, Ez, info) = MESTI(syst, B) returns the spatial field profiles
of Ex(x,y,z), Ey(x,y,z), and Ez(x,y,z) satisfying
[(∇ × ∇ ×) - (omega/c)^2*epsilon(x,y,z)]*[Ex(x,y,z); Ey(x,y,z); Ez(x,y,z)] =
i*omega*mu_0*[Jx(x,y,z); Jy(x,y,z); Jz(x,y,z)], where
The relative permittivity profile epsilon(x,y,z), frequency omega, and boundary conditions
are specified by structure "syst".
Note that relative permittivity profile epsilon(x,y,z) is a rank 2 tenor:
[epsilon_xx, epsilon_xy, epsilon_xz;
epsilon_yx, epsilon_yy, epsilon_yz;
epsilon_zx, epsilon_zy, epsilon_zz] in general.
Users can specify the diagonal terms only (epsilon_xx(x,y,z), epsilon_yy(x,y,z), and epsilon_zz(x,y,z))
or all of them.
Each column of matrix "B" specifies a distinct input source profile.
The returned "Ex", "Ey", "Ez" is a 4D array, such as Ex(:,:,:,i),
being the field profile Ex given the i-th input source profile. Same data structure for Ey and Ez.
The information of the computation is returned in structure "info".
MESTI uses finite-difference discretization on the Yee lattice, after which
the differential operator becomes an (nt_Ex*nt_Ey*nt_Ez)-by-(nt_Ex*nt_Ey*nt_Ez) sparse matrix A
where nt_Ex is the number of total grid sites for Ex and nt_Ex = nx_Ex*ny_Ex*nz_Ex
(nx_Ex, ny_Ex, nz_Ex) is the number of sites Ex is discretized onto. Same notation for Ey and Ez.
Ex = reshape((inv(A)*B)[nt_Ex,:], nx_Ex, ny_Ex, nz_Ex, :).
Ey = reshape((inv(A)*B)[nt_Ey,:], nx_Ey, ny_Ey, nz_Ey, :).
Ez = reshape((inv(A)*B)[nt_Ez,:], nx_Ez, ny_Ez, nz_Ez, :).
---2D TM field profile---
(Ex, info) = MESTI(syst, B) returns the spatial field profiles
of Ex(y,z) for 2D transverse-magnetic (TM) fields satisfying
[- (d/dy)^2 - (d/dz)^2 - (omega/c)^2*epsilon(y,z)] Ex(y,z) = i*omega*mu_0*Jx(y,z).
The returned 'Ex' is a 3D array, with Ex(:,:,i) being the field profile
of Ex given the i-th input source profile. The information of the computation is returned in structure 'info'.
MESTI uses finite-difference discretization on the Yee lattice, after which
the differential operator becomes an (ny_Ex*nz_Ex)-by-(ny_Ex*nz_Ex) sparse matrix A
where (ny_Ex, nz_Ex) is the number of sites Ex is discretized onto, and
Ex = reshape(inv(A)*B, ny_Ex, nz_Ex, []).
---Generalized scattering matrix S---
(S, info) = MESTI(syst, B, C) returns S = C*inv(A)*B where the solution
inv(A)*B is projected onto the output channels or locations of interest
through matrix C; each row of matrix "C" is a distinct output projection
profile, discretized into a 1-by-(nt_Ex*nt_Ey*nt_Ez) vector for 3D or
1-by-(ny_Ex*nz_Ex) vector for 2D TM fields in the same order as matrix A.
When the MUMPS is available, this is done by computing the Schur complement of an
augmented matrix K = [A,B;C,0] through a partial factorization.
(S, info) = MESTI(syst, B, C, D) returns S = C*inv(A)*B - D. This can be
used for the computation of scattering matrices, where S is the scattering
matrix, and matrix D can be derived analytically or computed as D =
C*inv(A0)*B - S0 from a reference system A0 for which the scattering matrix
S0 is known.
(Ex, Ey, Ez, info) = MESTI(syst, B, opts),
(Ex, info) = MESTI(syst, B, opts),
(S, info) = MESTI(syst, B, C, opts), and
(S, info) = MESTI(syst, B, C, D, opts) allow detailed options to be
specified with structure "opts".
MESTI only considers nonmagnetic materials.
This file checks and parses the parameters, and it can build matrices B and
C from its nonzero elements specified by the user (see details below). It
calls function mesti_build_fdfd_matrix() to build matrix A and function
mesti_matrix_solver!() to compute C*inv(A)*B or inv(A)*B, where most of the
computation is done.
=== Input Arguments ===
syst (Syst struct; required):
A structure that specifies the system, used to build the FDFD matrix A.
It contains the following fields:
syst.epsilon_xx (numeric array or matrix, real or complex, required):
For 3D systems, an nx_Ex-by-ny_Ex-by-nz_Ex array discretizing the relative permittivity
profile epsilon_xx(x,y,z). Specifically, syst.epsilon_xx(n,m,l) is the scalar
epsilon_xx(n,m,l) averaged over a cube with volume (syst.dx)^3 centered at
the point (x_n, y_m, z_l) where Ex(x,y,z) is located on the Yee lattice.
For 2D TM fields, an ny_Ex-by-nz_Ex matrix discretizing the relative permittivity
profile epsilon_xx(y,z).
syst.epsilon_xy (numeric array or matrix, real or complex, optional):
For 3D, an nx_Ez-by-ny_Ez-by-nz_Ex array discretizing the relative permittivity
profile epsilon_xy(x,y,z). Specifically, syst.epsilon_xy(n,m,l) is the scalar
epsilon_xy(n,m,l) averaged over a cube with volume (syst.dx)^3 centered at
the low corner on the Yee lattice (n,m,l).
syst.epsilon_xz (numeric array or nothing, real or complex, optional):
Discretizing the relative permittivity profile epsilon_xz(x,y,z),
analogous to syst.epsilon_xy.
syst.epsilon_yx (numeric array or nothing, real or complex, optional):
Discretizing the relative permittivity profile epsilon_yx(x,y,z),
analogous to syst.epsilon_xy.
syst.epsilon_yy (numeric array or nothing, real or complex, required for 3D):
Discretizing the relative permittivity profile epsilon_yy(x,y,z),
analogous to syst.epsilon_xx.
syst.epsilon_yz (numeric array or nothing, real or complex, optional):
Discretizing the relative permittivity profile epsilon_yz(x,y,z),
analogous to syst.epsilon_xy.
syst.epsilon_zx (numeric array or nothing, real or complex, optional):
Discretizing the relative permittivity profile epsilon_zx(x,y,z),
analogous to syst.epsilon_xy.
syst.epsilon_zy (numeric array or nothing, real or complex, optional):
Discretizing the relative permittivity profile epsilon_zy(x,y,z),
analogous to syst.epsilon_xy.
syst.epsilon_zz (numeric array or nothing, real or complex, required for 3D):
Discretizing the relative permittivity profile epsilon_zz(x,y,z),
analogous to syst.epsilon_xx.
syst.length_unit (string; optional):
Length unit, such as micron, nm, or some reference wavelength. This
code only uses dimensionless quantities, so syst.length_unit is never
used. This syst.length_unit is meant to help the user interpret the
units of (x,y,z), dx, wavelength, kx_B, ky_B, kz_B, etc.
syst.wavelength (numeric scalar, real or complex; required):
Vacuum wavelength 2*pi*c/omega, in units of syst.length_unit.
syst.dx (positive scalar; required):
Discretization grid size, in units of syst.length_unit.
syst.PML (PML structure or a vector of PML structure; optional):
Parameters of the perfectly matched layer (PML) used to simulate an
open boundary. Note that PML is not a boundary condition; it is a
layer placed within the simulation domain (just before the boundary)
that attenuates outgoing waves with minimal reflection.
In mesti(), the PML starts from the interior of the system
specified by syst.epsilon_ij (i = x,y,z and j = x,y,z),
and ends at the first or last pixel inside syst.epsilon_ij
(i = x,y,z and j = x,y,z). (Note: this is different from the function
mesti2s() that handles two-sided geometries, where the homogeneous spaces
on the low and high are specified separately through syst.epsilon_low
and syst.epsilon_high, and where PML is placed in such homogeneous space,
outside of the syst.epsilon_ij (i = x,y,z and j = x,y,z).
When only one set of PML parameters is used in the system (as is
the most common), such parameters can be specified with a scalar PML
structure syst.PML that contains the following fields:
npixels (positive integer scalar; required):
Number of PML pixels.
Note this is within syst.epsilon_ij (i = x,y,z and j = x,y,z),
not in addition to.
direction (string; required):
Direction(s) where PML is placed. Available choices are (case-insensitive):
"all" - (default) PML in x, y, and z directions for 3D and PML in y and z directions for 2D TM
"x" - PML in x direction
"y" - PML in y direction
"z" - PML in z direction
"xy" - PML in x and y directions
"xz" - PML in x and z directions
"yz" - PML in y and z directions
"yx" - PML in x and y directions
"zx" - PML in x and z directions
"zy" - PML in y and z directions
side (string; optional):
Side(s) where PML is placed.Available choices are (case-insensitive):
"both" - (default) PML on both sides
"-" - one-sided PML; end at the first pixel
"+" - one-sided PML; end at the last pixel
power_sigma (non-negative scalar; optional):
Power of the polynomial grading for the conductivity sigma;
defaults to 3.
sigma_max_over_omega (non-negative scalar; optional):
Conductivity at the end of the PML; By default, it is
set to an optimized value based on resolution and the
background refractive index. This is used to attenuate propagating waves.
power_kappa (non-negative scalar; optional):
Power of the polynomial grading for the real-coordinate-stretching
factor kappa; defaults to 3.
kappa_max (real scalar no smaller than 1; optional):
Real-coordinate-stretching factor at the end of the PML; By default,
it is set to an optimized value based on resolution and the
background refractive index. This is used to accelerate the attenuation
of evanescent waves. kappa_max = 1 means no real-coordinate stretching.
power_alpha (non-negative scalar; optional):
Power of the polynomial grading for the CFS alpha factor;
defaults to 1.
alpha_max_over_omega (non-negative scalar; optional):
Complex-frequency-shifting (CFS) factor at the beginning
of the PML. This is typically used in time-domain simulations
to suppress late-time (low-frequency) reflections.
We don't use it by default (alpha_max_over_omega = 0)
since we are in frequency domain.
We use the following PML coordinate-stretching factor:
s(p) = kappa(p) + sigma(p)./(alpha(p) - i*omega)
with
sigma(p)/omega = sigma_max_over_omega*(p.^power_sigma),
kappa(p) = 1 + (kappa_max-1)*(p.^power_kappa),
alpha(p)/omega = alpha_max_over_omega*((1-p).^power_alpha),
where omega is frequency, and p goes linearly from 0 at the beginning
of the PML to 1 at the end of the PML.
If isdefined(syst, :PML) = false, which means no PML on any side. PML is
only placed on the side(s) specified by syst.PML.
When multiple sets of PML parameters are used in the system (e.g.,
a thinner PML on one side, a thicker PML on another side), these
parameters can be specified with a vector of PML strcture.
syst.PML = [PML_1, PML_2, ...],
with PML_1 and PML_2 each being a structure containing the above
fields; they can specify different PML parameters on different sides.
Each side cannot be specified more than once.
With real-coordinate stretching, PML can attenuate evanescent waves
more efficiently than free space, so there is no need to place free
space in front of PML.
The PML thickness should be chosen based on the acceptable level of
reflectivity given the discretization resolution and the range of wave
numbers (i.e., angles) involved; more PML pixels gives lower
reflectivity. Typically 10-40 pixels are sufficient.
syst.PML_type (string; optional):
Type of PML. Available choices are (case-insensitive):
"UPML" - (default) uniaxial PML
"SC-PML" - stretched-coordinate PML
The two are mathematically equivalent, but using SC-PML has lower
condition number.
syst.xBC (string or nothing; optional):
Boundary condition (BC) at the two ends in x direction, effectively
specifying Ex(n,m,l) at n=0 and n=nx_Ex+1,
Ey(n,m,l) at n=0 and n=nx_Ey+1,
Ez(n,m,l) at n=0 and n=nx_Ez+1.
one pixel beyond the computation domain. Available choices are:
"Bloch" - Ex(n+nx_Ex,m,l) = Ex(n,m,l)*exp(1i*syst.kx_B*nx_Ex*syst.dx),
Ey(n+nx_Ey,m,l) = Ey(n,m,l)*exp(1i*syst.kx_B*nx_Ey*syst.dx),
Ez(n+nx_Ez,m,l) = Ez(n,m,l)*exp(1i*syst.kx_B*nx_Ez*syst.dx).
"periodic" - equivalent to "Bloch" with syst.kx_B = 0
"PEC" - Ex(0,m,l) = Ex(1,m,l); Ex(nx_Ex+1,m,l) = Ez(nx_Ex,m,l),
Ey(0,m,l) = Ey(nx_Ey+1,m,l) = 0,
Ez(0,m,l) = Ez(nx_Ez+1,m,l) = 0.
"PMC" - Ex(0,m,l) = Ex(nx_Ex+1,m,l) = 0,
Ey(0,m,l) = Ey(1,m,l); Ey(nx_Ey+1,m,l) = Ey(nx_Ey,m,l),
Ez(0,m,l) = Ez(1,m,l); Ez(nx_Ez+1,m,l) = Ez(nx_Ez,m,l).
"PECPMC" - Ex(0,m,l) = Ex(1,m,l); Ex(nx_Ex+1,m,l) = 0,
Ey(0,m,l) = 0; Ey(nx_Ey+1,m,l) = Ey(nx_Ey,m,l),
Ez(0,m,l) = 0; Ez(nx_Ez+1,m,l) = Ez(nx_Ez,m,l),
"PMCPEC" - Ex(0,m,l) = 0; Ex(nx_Ex+1,m,l) = Ex(nx_Ex,m,l),
Ey(0,m,l) = Ey(1,m,l); Ey(nx_Ey+1,m,l) = 0,
Ez(0,m,l) = Ez(1,m,l); Ez(nx_Ez+1,m,l) = 0.
where PEC stands for perfect electric conductor and PMC stands for perfect
magnetic conductor.
By default, syst.xBC = "Bloch" if syst.kx_B is given; otherwise syst.xBC = "PEC"
The choice of syst.xBC has little effect on the numerical accuracy
when PML is used.
For 2D TM fields, xBC = nothing.
syst.yBC (string; optional):
Boundary condition in y direction, analogous to syst.xBC.
syst.zBC (string; optional):
Boundary condition in z direction, analogous to syst.xBC.
syst.kx_B (real scalar or nothing; optional):
Bloch wave number in x direction, in units of 1/syst.length_unit.
syst.kx_B is only used when syst.xBC = "Bloch". It is allowed to
specify a complex-valued syst.kx_B, but a warning will be displayed.
For 2D TM fields, kx_B = nothing.
syst.ky_B (real scalar; optional):
Bloch wave number in y direction, analogous to syst.kx_B.
syst.kz_B (real scalar; optional):
Bloch wave number in z direction, analogous to syst.kx_B.
B (numeric matrix or vector of source structure; required):
Matrix specifying the input source profiles B in the C*inv(A)*B - D or
C*inv(A)*B or inv(A)*B returned. When the input argument B is a matrix,
it is directly used, and size(B,1) must equal nt_Ex*nt_Ey*nt_Ez for 3D system and nt_Ex for 2D TM case;
each column of B specifies a source profile, placed on the grid points of E.
Note that matrix A is (syst.dx)^2 times the differential operator and
is unitless, so each column of B is (syst.dx)^2 times the i*omega*mu_0*[Jx;Jy;Jz] on
the right-hand side of the differential equation and has the same unit as E.
Instead of specifying matrix B directly, one can specify only its
nonzero parts, from which mesti() will build the sparse matrix B. To do
so, B in the input argument should be set as a source structure vector;
(Bx_struct, By_struct, Bz_struct) here we refer to such source structure
as Bx_struct as source for x component to distinguish it from the
resulting matrix B. If for every column of matrix B, all of its nonzero
elements are spatially located within a cuboid (e.g., block sources),
one can use the following fields:
Bx_struct.pos (six-element integer vector or four-element integer vector): Bx_struct.pos =
[n1, m1, l1, d, w, h] specifies the location and the size of the
cuboid on Ex-grid. Here, (n1, m1, l1) is the index of the (x, y, z) coordinate of
the smaller-x, smaller-y, and smaller-z corner of the cuboid,
at the location of p(n1, m1, l1) where p = Ex; (d, w, h) is the depth, width, and height
of the cuboid, such that (n2, m2, l2) = (n1+d-1, m1+w-1, l1+h-1) is the index
of the higher-index corner of the cuboid.
For 2D TM case, users only need to specify [m1, l1, w, h].
Bx_struct.data (2D, 3D or 4D numeric array): nonzero elements of matrix B for x component
within the cuboid specified by Bx_struct.pos.
When it is a 4D array (general 3D system), Bx_struct.data[n',m',l',a] is the a-th input
source at the location of p(n=n1+n'-1, m=m1+m'-1, l=l1+l'-1), which becomes
B[n+(m-1)*nx_Ex+(l-1)*nx_Ex*ny_Ex, a]. In other words, Bx_struct.data[:,:,:,a] gives the
sources at the cuboid p(n1+(0:(d-1)), m1+(0:(w-1)), l1+(0:(h-1))). So,
size(Bx_struct.data) must equal (d, w, h, number of inputs).
When it is a 3D array (2D TM system), Bx_struct.data(m',l',a) is the a-th input
source at the location of p(m=m1+m'-1, l=l1+l'-1), which becomes
Bx(m+(l-1)*ny, a). In other words, Bx_struct.data(:,:,a) gives the
sources at the rectangle p(m1+(0:(w-1)), l1+(0:(h-1))). So,
size(Bx_struct.data, [1,2]) must equal [w, h], and
size(Bx_struct.data, 3) is the number of inputs.
Alternatively, Bx_struct.data can be a 2D array that is
equivalent to reshape(data_in_4D_array, d*w*h, :) or reshape(data_in_3D_array, w*h, :), in which case
size(Bx_struct.data, 2) is the number of inputs; in this case,
Bx_struct.data can be a sparse matrix, and its sparsity will be
preserved when building matrix B.
By_struct.pos (six-element integer vector): By_struct.pos =
[n1, m1, l1, d, w, h] specifies the location and the size of the
cuboid on Ey-grid, analogous to Bx_struct.pos.
By_struct.data (2D or 4D numeric array): nonzero elements of matrix B for y component
within the cuboid specified by By_struct.pos, analogous to Bx_struct.data.
Bz_struct.pos (six-element integer vector): Bz_struct.pos =
[n1, m1, l1, d, w, h] specifies the location and the size of the
cuboid on Ez-grid, analogous to Bx_struct.pos.
Bz_struct.data (2D or 4D numeric array): nonzero elements of matrix B for z component
within the cuboid specified by Bz_struct.pos, analogous to Bx_struct.data.
If different inputs are located within different cuboids (e.g.,
inputs from plane sources on the low and separate inputs from plane
sources on the high), Bx_struct can be a structure array with multiple
elements (e.g., Bx_struct.pos[1] and Bx_struct.data[1] specify plane sources
on the low; Bx_struct.pos[2] and Bx_struct.data[2] specify plane sources on
the high); these inputs are treated separately, and the total number of
inputs is size(Bx_struct.data[1], 4) + size(Bx_struct.data[2], 4) + ... +
size(Bx_struct.data[end], 4).
If the nonzero elements of matrix B do not have cuboid shapes in
space [e.g., for total-field/scattered-field (TF/SF) simulations], one
can use a source structure with the following fields:
Bx_struct.ind (integer vector): linear indices of the spatial
locations of the nonzero elements of matrix B for x component, such that
p(Bx_struct.ind) are the points where the source is placed. Such linear
indices can be constructed from Base._sub2ind().
Bx_struct.data (2D numeric matrix): nonzero elements of matrix B for x component
at the locations specified by Bx_struct.ind. Specifically,
Bx_struct.data[i,a] is the a-th input source at the location of
p(Bx_struct.ind[i]), which becomes B[Bx_struct.ind[i], a]. So,
size(Bx_struct.data, 1) must equal length(Bx_struct.ind), and
size(Bx_struct.data, 2) is the number of inputs.
By_struct.ind (integer vector): linear indices of the spatial
locations of the nonzero elements of matrix B for y component, analogous to Bx_struct.ind.
By_struct.data (2D numeric matrix): nonzero elements of matrix B for y component
at the locations specified by By_struct.ind, analogous to Bx_struct.data.
Bz_struct.ind (integer vector): linear indices of the spatial
locations of the nonzero elements of matrix B for z component, analogous to Bx_struct.ind.
Bz_struct.data (2D numeric matrix): nonzero elements of matrix B for z component
at the locations specified by Bz_struct.ind, analogous to Bx_struct.data.
Similarly, one can use Bx_struct.ind[1], Bx_struct.ind[2] etc together
with Bx_struct.data[1], Bx_struct.data[2] etc to specify inputs at
different sets of locations. Every element of the structure array must
have the same fields (e.g., one cannot specify Bx_struct.pos[1] and
Bx_struct.ind[2]), so the more general Bx_struct.ind syntax should be used
when some of the inputs are rectangular and some are not.
C (numeric matrix, vector of source structure, or "transpose(B)"; optional):
Matrix specifying the output projections in the C*inv(A)*B - D or
C*inv(A)*B returned. When the input argument C is a matrix, it is
directly used, and size(C,2) must equal nt_Ex*nt_Ey*nt_Ez for 3D system and nt_Ex for 2D TM case;
each row of C specifies a projection profile, placed on the grid points of E.
Scattering matrix computations often have C = transpose(B); if that
is the case, the user can set C = "transpose(B)" as a string,
and it will be replaced by transpose(B) in the code.
For field-profile computations, the user can simply omit C from the
input arguments, as in mesti(syst, B). If opts is needed, the user can use
mesti(syst, B, opts).
Similar to B, here one can specify only the nonzero parts of the
output matrix C, from which mesti() will build the sparse matrix C. The
syntax of vector of source structure (Cx_struct, Cy_struct, Cz_struct) is
the same as for B, summarized below. If for every row of matrix
C, all of its nonzero elements are spatially located withing a cuboid
(e.g., projection of fields on a line), one can set the input argument C
to be a structure array (referred to as Cx_struct below) with the
following fields:
Cx_struct.pos (six-element integer vector or four-element integer vector): Cx_struct.pos =
[n1, m1, l1, d, w, h] specifies the location and the size of the
cuboid on Ex-grid. Here, (n1, m1, l1) is the index of the (x, y, z) coordinate of
the smaller-x, smaller-y, and smaller-z corner of the cuboid,
at the location of p(n1, m1, l1) where p = Ex; (d, w, h) is the depth, width, and height
of the cuboid, such that (n2, m2, l2) = (n1+d-1, m1+w-1, l1+h-1) is the index
of the higher-index corner of the cuboid.
For 2D TM case, users only need to specify [m1, l1, w, h].
Cx_struct.data (2D, 3D or 4D numeric array): nonzero elements of matrix C for x component
within the cuboid specified by Cx_struct.pos.
When it is a 4D array (general 3D system), Cx_struct.data[n',m',l',a] is the a-th input
source at the location of p(n=n1+n'-1, m=m1+m'-1, l=l1+l'-1), which becomes
C[n+(m-1)*nx_Ex+(l-1)*nx_Ex*ny_Ex, a]. In other words, Cx_struct.data[:,:,:,a] gives the
sources at the cuboid p(n1+(0:(d-1)), m1+(0:(w-1)), l1+(0:(h-1))). So,
size(Cx_struct.data) must equal (d, w, h, number of inputs).
When it is a 3D array (2D system), Cx_struct.data(m',l',a) is the a-th input
source at the location of p(m=m1+m'-1, l=l1+l'-1), which becomes
Cx(m+(l-1)*ny, a). In other words, Cx_struct.data(:,:,a) gives the
sources at the rectangle p(m1+(0:(w-1)), l1+(0:(h-1))). So,
size(Cx_struct.data, [1,2]) must equal [w, h], and
size(Cx_struct.data, 3) is the number of inputs.
Alternatively, Cx_struct.data can be a 2D array that is
equivalent to reshape(data_in_4D_array, d*w*h, :) or reshape(data_in_3D_array, w*h, :), in which case
size(Cx_struct.data, 2) is the number of inputs; in this case,
Cx_struct.data can be a sparse matrix, and its sparsity will be
preserved when building matrix C.
Cy_struct.pos (six-element integer vector): Cy_struct.pos =
[n1, m1, l1, d, w, h] specifies the location and the size of the
cuboid on Ey-grid, analogous to Cx_struct.pos.
Cy_struct.data (2D or 4D numeric array): nonzero elements of matrix C for y component
within the cuboid specified by Cy_struct.pos, analogous to Cx_struct.data.
Cz_struct.pos (six-element integer vector): Cz_struct.pos =
[n1, m1, l1, d, w, h] specifies the location and the size of the
cuboid on Ez-grid, analogous to Cx_struct.pos.
Cz_struct.data (2D or 4D numeric array): nonzero elements of matrix C for z component
within the cuboid specified by Cz_struct.pos, analogous to Cx_struct.data.
If the nonzero elements of matrix C do not have rectangular shapes in
space [e.g., for near-field-to-far-field transformations], one can set C
to a structure array with the following fields:
Cx_struct.ind (integer vector): linear indices of the spatial
locations of the nonzero elements of matrix C for x component, such that
p(Cx_struct.ind) are the points where the source is placed.
Cx_struct.data (2D numeric matrix): nonzero elements of matrix C for x component
at the locations specified by Cx_struct.ind. Specifically,
Cx_struct.data[i,a] is the a-th input source at the location of
p(Cx_struct.ind[i]), which becomes C[Cx_struct.ind[i], a]. So,
size(Cx_struct.data, 1) must equal length(Cx_struct.ind), and
size(Cx_struct.data, 2) is the number of inputs.
Cy_struct.ind (integer vector): linear indices of the spatial
locations of the nonzero elements of matrix C for y component, analogous to Cx_struct.ind.
Cy_struct.data (2D numeric matrix): nonzero elements of matrix C for y component
at the locations specified by Cy_struct.ind, analogous to Cx_struct.data.
Cz_struct.ind (integer vector): linear indices of the spatial
locations of the nonzero elements of matrix C for z component, analogous to Cx_struct.ind.
Cz_struct.data (2D numeric matrix): nonzero elements of matrix C for z component
at the locations specified by Cz_struct.ind, analogous to Cx_struct.data.
Like in Bx_struct, one can use structure arrays with multiple elements
to specify outputs at different spatial locations.
D (numeric matrix; optional):
Matrix D in the C*inv(A)*B - D returned, which specifies the baseline
contribution; size(D,1) must equal size(C,1), and size(D,2) must equal
size(B,2). When D is not an input argument, it will not be subtracted from C*inv(A)*B.
For field-profile computations where C is not an input argument, D should not be an input
argument, either.
opts (Opts structure; optional):
A structure that specifies the options of computation.
It can contain the following fields (all optional):
opts.verbal (boolean scalar; optional, defaults to true):
Whether to print system information and timing to the standard output.
opts.prefactor (numeric scalar, real or complex; optional):
When opts.prefactor is given, mesti() will return
opts.prefactor*C*inv(A)*B - D or opts.prefactor*C*inv(A)*B or
opts.prefactor*inv(A)*B. Such prefactor makes it easier to use C =
transpose(B) to take advantage of reciprocity. Defaults to 1.
opts.exclude_PML_in_field_profiles (boolean scalar; optional, defaults to false):
When opts.exclude_PML_in_field_profiles = true, the PML pixels
(specified by syst.PML.npixels) are excluded from the returned
field_profiles on each side where PML is used; otherwise the full
field profiles are returned. Only used for field-profile computations
(i.e., when the output projection matrix C is not given).
opts.solver (string; optional):
The solver used for sparse matrix factorization. Available choices
are (case-insensitive):
"MUMPS" - (default when MUMPS is available) To use MUMPS, MUMPS must
be installed and the Julia environment variable "MUMPS_PREFIX"
should be specified.
"JULIA" - (default when MUMPS is not available) Uses the built-in
lu() function in JULIA, which uses UMFPACK.
MUMPS is faster and uses less memory than lu(), and is required for
the APF method.
opts.method (string; optional):
The solution method. Available choices are (case-insensitive):
"APF" - Augmented partial factorization. C*inv(A)*B is obtained
through the Schur complement of an augmented matrix
K = [A,B;C,0] using a partial factorization. Must have
opts.solver = "MUMPS". This is the most efficient method,
but it cannot be used for computing the full field profile
X=inv(A)*B or with iterative refinement.
"FG" - Factorize and group. Factorize A=L*U, and obtain C*inv(A)*B
through C*inv(U)*inv(L)*B with optimized grouping. Must
have opts.solver = "JULIA". This is slightly better than
"FS" when MUMPS is not available, but it cannot be used for
computing the full field profile X=inv(A)*B.
"FS" - Factorize and solve. Factorize A=L*U, solve for X=inv(A)*B
with forward and backward substitutions, and project with
C as C*inv(A)*B = C*X. Here, opts.solver can be either
"MUMPS" or "JULIA", and it can be used for computing
the full field profile X=inv(A)*B or with iterative
refinement.
"C*inv(U)*inv(L)*B" - Same as "FG".
"factorize_and_solve" - Same as "FS".
By default, if C is given and opts.iterative_refinement = false, then
"APF" is used when opts.solver = "MUMPS", and "C*inv(U)*inv(L)*B" is
used when opts.solver = "JULIA". Otherwise, "factorize_and_solve" is
used.
opts.clear_BC (boolean scalar; optional, defaults to false):
When opts.clear_BC = true, variables "B" and "C" will be cleared in
the caller's workspace to reduce peak memory usage. Can be used when B
and/or C take up significant memory and are not needed after calling
mesti(). However, it is not implemented and does not work for current
julia version.
opts.clear_syst (boolean scalar; optional, defaults to false):
When opts.clear_syst = true, variable "syst" will be cleared in the
caller's workspace to reduce peak memory usage. This can be used when
syst.epsilon_xx, syst.epsilon_yy, and syst.epsilon_zz take up
significant memory and are not needed after calling mesti(). However,
it is not implemented and does not work for current julia version.
opts.clear_memory (boolean scalar; optional, defaults to true):
Whether or not to clear variables inside mesti() to reduce peak memory
usage. If opts.clear_memory == true, we call GC.gc() inside mesti().
opts.verbal_solver (boolean scalar; optional, defaults to false):
Whether to have the solver print detailed information to the standard
output. Note the behavior of output from MUMPS depends on compiler.
opts.use_single_precision_MUMPS (boolean scalar; optional, defaults to true):
Whether to use single precision version of MUMPS; used only when
opts.solver = "MUMPS". Using single precision version of MUMPS can
reduce memory usage and computing time.
opts.use_METIS (boolean scalar; optional, defaults to false in 2D and to true in 3D):
Whether to use METIS (instead of the default AMD) to compute the
ordering in MUMPS. Using METIS can sometimes reduce memory usage
and/or factorization and solve time, but it typically takes longer at
the analysis (i.e., ordering) stage in 2D. In 3D METIS is general better
than AMD.
opts.nrhs (positive integer scalar; optional):
The number of right-hand sides (number of columns of the input matrix
B) to consider simultaneously, used only when opts.method =
"factorize_and_solve" and C is given. Defaults to 1 if
opts.iterative_refinement = true, 10 if opts.solver = "MUMPS" with
opts.iterative_refinement = false, 4 otherwise.
opts.store_ordering (boolean scalar; optional, defaults to false):
Whether to store the ordering sequence (permutation) for matrix A or
matrix K; only possible when opts.solver = "MUMPS". If
opts.store_ordering = true, the ordering will be returned in
info.ordering.
opts.ordering (positive integer vector; optional):
A user-specified ordering sequence for matrix A or matrix K, used only
when opts.solver = "MUMPS". Using the ordering from a previous
computation can speed up (but does not eliminate) the analysis stage.
The matrix size must be the same, and the sparsity structure should be
similar among the previous and the current computation.
opts.analysis_only (boolean scalar; optional, defaults to false):
When opts.analysis_only = true, the factorization and solution steps
will be skipped, and S = nothing will be returned. The user can use
opts.analysis_only = true with opts.store_ordering = true to return
the ordering for A or K; only possible when opts.solver = "MUMPS".
opts.nthreads_OMP (positive integer scalar; optional):
Number of OpenMP threads used in MUMPS; overwrites the OMP_NUM_THREADS
environment variable.
opts.use_L0_threads (logical scalar; optional,
defaults to true in 1D/2D/2.5D and to false in full 3D):
If MUMPS is multithread, whether to use tree parallelism (so-called
L0-threads layer) in MUMPS. Please refer to Sec. 5.23 'Improved
multithreading using tree parallelism' in MUMPS 5.7.1 Users' guide.
This typically enhances the time performance, but marginally increases
the memory usage. In full 3D ((width in x)*(width in y)/(thickness in z) < 100),
memory usage is critical and we set it false. Otherwise, it is enabled
by default.
opts.write_LU_factor_to_disk (logical scalar; optional, defaults to true):
An out-of-core (disk is used as an extension to main memory) facility
is utilized to write the complete matrix of factors to disk in the
factorization phase and read them in the solve phase. This can significantly
reduce the memory requirement while not increasing much the factorization time.
The extra cost of the out-of-core feature is thus mainly during the solve phase,
where factors have to be read from disk.
opts.iterative_refinement (boolean scalar; optional, defaults to false):
Whether to use iterative refinement in MUMPS to lower round-off
errors. Iterative refinement can only be used when opts.solver =
"MUMPS" and opts.method = "factorize_and_solve" and C is given, in
case opts.nrhs must equal 1. When iterative refinement is used, the
relevant information will be returned in info.itr_ref_nsteps,
info.itr_ref_omega_1, and info.itr_ref_omega_2.
opts.use_BLR (logical scalar; optional, defaults to false):
Whether to use block low-rank approximation in MUMPS to possibly lower computational
cost (but in most case it does not). It can only be used when opts.solver = "MUMPS".
opts.threshold_BLR (positive real scalar; optional):
The dropping parameter controls the accuracy of the block low-rank approximations.
It can only be used when opts.solver = "MUMPS" and opts.use_BLR = true.
Please refer to the section of BLR API in MUMPS userguide.
opts.icntl_36 (positive integer scalar; optional):
It controls the choice of the BLR factorization variant.
It can only be used when opts.solver = "MUMPS" and opts.use_BLR = true.
Please refer to the section of BLR API in MUMPS userguide.
opts.icntl_38 (positive integer scalar; optional):
It estimated compression rate of LU factors.
It can only be used when opts.solver = "MUMPS" and opts.use_BLR = true.
Please refer to the section of BLR API in MUMPS userguide.
=== Output Arguments ===
---When users specify C---
S (full numeric matrix):
The generalized scattering matrix S = C*inv(A)*B or S = C*inv(A)*B - D.
---or when users do not specify C---
When opts.exclude_PML_in_field_profiles = false,
For 3D system,
Ex (4D array):
Electrical field profile for Ex component
(nx_Ex, ny_Ex, nz_Ex, number of inputs) = size(Ex)
Ey (4D array):
Electrical field profile for Ey component
(nx_Ey, ny_Ey, nz_Ey, number of inputs) = size(Ey)
Ez (4D array):
Electrical field profile for Ez component
(nx_Ez, ny_Ez, nz_Ez, number of inputs) = size(Ez)
---For 2D TM case---
Ex (3D array):
Electrical field profile for Ex component
(ny_Ex, nz_Ex, number of inputs) = size(Ex)
When opts.exclude_PML_in_field_profiles = true, the PML pixels
(specified by syst.PML.npixels) are excluded from Ex, Ey, and Ez on each
side where PML is used.
info (Info structure):
A structure that contains the following fields:
info.opts (Opts structure):
The final "opts" used, excluding the user-specified matrix ordering.
info.timing_init (non-negative scalar):
Timing in the "init" stages
info.timing_build (non-negative scalar):
Timing in the "building" stages
info.timing_analyze (non-negative scalar):
Timing in the "analysis" stages
info.timing_factorize (non-negative scalar):
Timing in the "factorizing" stages
info.timing_total (non-negative scalar):
Timing in total
info.xPML (two-element cell array; optional);
PML parameters on the sides of x direction, if used.
info.yPML (two-element cell array; optional);
PML parameters on the sides of y direction, if used.
info.zPML (two-element cell array; optional);
PML parameters on the sides of z direction, if used.
info.ordering_method (string; optional):
Ordering method used in MUMPS.
info.ordering (positive integer vector; optional):
Ordering sequence returned by MUMPS when opts.store_ordering = true.
info.itr_ref_nsteps (integer vector; optional):
Number of steps of iterative refinement for each input, if
opts.iterative_refinement = true; 0 means no iterative refinement.
info.itr_ref_omega_1 (real vector; optional):
Scaled residual omega_1 at the end of iterative refinement for each
input; see MUMPS user guide section 3.3.2 for definition.
info.itr_ref_omega_2 (real vector; optional):
Scaled residual omega_2 at the end of iterative refinement for each
input; see MUMPS user guide section 3.3.2 for definition.
See also: mesti_build_fdfd_matrix, mesti_matrix_solver!, mesti2s
"""
function mesti(syst::Syst, B::Union{SparseMatrixCSC{Int64,Int64},SparseMatrixCSC{Float64, Int64},SparseMatrixCSC{ComplexF64,Int64},Array{Int64,2},Array{Float64,2},Array{ComplexF64,2},Vector{Source_struct}}, C::Union{SparseMatrixCSC{Int64,Int64},SparseMatrixCSC{Float64, Int64},SparseMatrixCSC{ComplexF64,Int64},Array{Int64,2},Array{Float64,2},Array{ComplexF64,2}, Vector{Source_struct},String,Nothing}, D::Union{SparseMatrixCSC{Int64,Int64},SparseMatrixCSC{Float64, Int64},SparseMatrixCSC{ComplexF64,Int64},Array{Int64,2},Array{Float64,2},Array{ComplexF64,2},Array{Int32,2},Array{Float32,2},Array{ComplexF32,2},Nothing}, opts::Union{Opts,Nothing})
if ~(stacktrace()[2].func == :mesti2s)
# Make deepcopy of them to avoid mutating input argument
syst = deepcopy(syst); opts = deepcopy(opts)
end
## Part 1.1: Check validity of syst, assign default values to its fields, and parse BC and PML specifications
t0 = time()
if ~isdefined(syst, :wavelength); throw(ArgumentError("Input argument syst must have field \"wavelength\".")); end
if ~isdefined(syst, :dx); throw(ArgumentError("Input argument syst must have field \"dx\".")); end
if ~(syst.dx > 0); throw(ArgumentError("syst.dx must be a positive scalar.")); end
# Take care of the 2D TM cases
if ndims(syst.epsilon_xx) == 2
use_2D_TM = true
if (isdefined(syst, :epsilon_yy) && ~isa(syst.epsilon_yy, Nothing)) || (isdefined(syst, :epsilon_zz) && ~isa(syst.epsilon_zz, Nothing)) || (isdefined(syst, :epsilon_xy) && ~isa(syst.epsilon_xy, Nothing)) || (isdefined(syst, :epsilon_xz) && ~isa(syst.epsilon_xz, Nothing)) || (isdefined(syst, :epsilon_yx) && ~isa(syst.epsilon_yx, Nothing)) || (isdefined(syst, :epsilon_yz) && ~isa(syst.epsilon_yz, Nothing)) || (isdefined(syst, :epsilon_zx) && ~isa(syst.epsilon_zx, Nothing)) || (isdefined(syst, :epsilon_zy) && ~isa(syst.epsilon_zy, Nothing))
throw(ArgumentError("Only syst.epsilon_xx is required for 2D TM fields Ex(y,z), but other components should not be given or should be nothing"))
end
else
use_2D_TM = false
end
if ~use_2D_TM
if (isdefined(syst, :epsilon_yy) && ~isa(syst.epsilon_yy, Nothing)) || (isdefined(syst, :epsilon_zz) && ~isa(syst.epsilon_zz, Nothing)) || (isdefined(syst, :epsilon_xy) && ~isa(syst.epsilon_xy, Nothing)) || (isdefined(syst, :epsilon_xz) && ~isa(syst.epsilon_xz, Nothing)) || (isdefined(syst, :epsilon_yx) && ~isa(syst.epsilon_yx, Nothing)) || (isdefined(syst, :epsilon_yz) && ~isa(syst.epsilon_yz, Nothing)) || (isdefined(syst, :epsilon_zx) && ~isa(syst.epsilon_zx, Nothing)) || (isdefined(syst, :epsilon_zy) && ~isa(syst.epsilon_zy, Nothing))
include_off_diagonal_epsilon = true
if ~isdefined(syst, :epsilon_xy); syst.epsilon_xy = nothing; end
if ~isdefined(syst, :epsilon_xz); syst.epsilon_xz = nothing; end
if ~isdefined(syst, :epsilon_yx); syst.epsilon_yx = nothing; end
if ~isdefined(syst, :epsilon_yz); syst.epsilon_yz = nothing; end
if ~isdefined(syst, :epsilon_zx); syst.epsilon_zx = nothing; end
if ~isdefined(syst, :epsilon_zy); syst.epsilon_zy = nothing; end
else
include_off_diagonal_epsilon = false
end
# Number of grid points in x, y, and z for Ex, Ey, and Ez
(nx_Ex, ny_Ex, nz_Ex) = size(syst.epsilon_xx)
(nx_Ey, ny_Ey, nz_Ey) = size(syst.epsilon_yy)
(nx_Ez, ny_Ez, nz_Ez) = size(syst.epsilon_zz)
# Total number of grid points for Ex, Ey, and Ez
nt_Ex = nx_Ex*ny_Ex*nz_Ex
nt_Ey = nx_Ey*ny_Ey*nz_Ey
nt_Ez = nx_Ez*ny_Ez*nz_Ez
else
# Consider 2D TM field: Ex(y,z)
# Number of grid points in y and z for Ex
(ny_Ex, nz_Ex) = size(syst.epsilon_xx)
# Total number of grid points for Ex
nt_Ex = ny_Ex*nz_Ex
end
# Check that the user did not accidentally use options only in mesti2s()
if isdefined(syst, :epsilon_low) && ~isa(syst.epsilon_low, Nothing)
@warn "syst.epsilon_low is not used in mesti(); will be ignored."
syst.epsilon_low = nothing
end
if isdefined(syst, :epsilon_high) && ~isa(syst.epsilon_high, Nothing)
@warn "syst.epsilon_high is not used in mesti(); will be ignored."
syst.epsilon_high = nothing
end
if isdefined(syst, :zPML) && ~isa(syst.zPML, Nothing)
@warn "syst.zPML is not used in mesti(); will be ignored."
syst.zPML = nothing
end
if ~use_2D_TM
# Check boundary condition in x
if isdefined(syst, :kx_B)
if isdefined(syst, :xBC) && lowercase(syst.xBC) != "bloch"
throw(ArgumentError("When syst.kx_B is given, syst.xBC must be \"Bloch\" if specified."))
end
syst.xBC = "Bloch"
# mesti_build_fdfd_matrix() uses (kx_B,ky_B,kz_B)*periodicity as the input arguments xBC, yBC, and zBC for Bloch BC
xBC = (syst.kx_B)*(nx_Ex*syst.dx) # dimensionless
else
# Defaults to Dirichlet boundary condition unless syst.kx_B is given
if ~isdefined(syst, :xBC)
syst.xBC = "PEC"
elseif ~(lowercase(syst.xBC) in ["bloch", "periodic", "pec", "pmc", "pecpmc", "pmcpec"])
throw(ArgumentError("syst.xBC = \"$(syst.xBC)\" is not a supported option; type ''? mesti'' for supported options."))
elseif lowercase(syst.xBC) == "bloch"
throw(ArgumentError("syst.xBC = \"Bloch\" but syst.kx_B is not given."))
end
xBC = syst.xBC
end
end
# Check boundary condition in y
if isdefined(syst, :ky_B)
if isdefined(syst, :yBC) && lowercase(syst.yBC) != "bloch"
throw(ArgumentError("When syst.ky_B is given, syst.yBC must be \"Bloch\" if specified."))
end
syst.yBC = "Bloch"
# mesti_build_fdfd_matrix() uses (kx_B,ky_B,kz_B)*periodicity as the input arguments xBC, yBC, and zBC for Bloch BC
if ~use_2D_TM
yBC = (syst.ky_B)*(ny_Ey*syst.dx) # dimensionless
else
yBC = (syst.ky_B)*(ny_Ex*syst.dx) # dimensionless
end
else
# Defaults to Dirichlet boundary condition unless syst.ky_B is given
if ~isdefined(syst, :yBC)
syst.yBC = "PEC"
elseif ~(lowercase(syst.yBC) in ["bloch", "periodic", "pec", "pmc", "pecpmc", "pmcpec"])
throw(ArgumentError("syst.yBC = \"$(syst.yBC)\" is not a supported option; type ''? mesti'' for supported options."))
elseif lowercase(syst.yBC) == "bloch"
throw(ArgumentError("syst.yBC = \"Bloch\" but syst.ky_B is not given."))
end
yBC = syst.yBC
end
# Check boundary condition in z
if isdefined(syst, :kz_B)
if isdefined(syst, :zBC) && lowercase(syst.zBC) != "bloch"
throw(ArgumentError("When syst.kz_B is given, syst.zBC must be \"Bloch\" if specified."))
end
syst.zBC = "Bloch"
# mesti_build_fdfd_matrix() uses (kx_B,ky_B,kz_B)*periodicity as the input arguments xBC, yBC, and zBC for Bloch BC
if ~use_2D_TM
zBC = (syst.kz_B)*(nz_Ez*syst.dx) # dimensionless
else
zBC = (syst.kz_B)*(nz_Ex*syst.dx) # dimensionless
end
else
# Defaults to Dirichlet boundary condition unless syst.kz_B is given
if ~isdefined(syst, :zBC)
syst.zBC = "PEC"
elseif ~(lowercase(syst.zBC) in ["bloch", "periodic", "pec", "pmc", "pecpmc", "pmcpec"])
throw(ArgumentError("syst.zBC = \"$(syst.zBC)\" is not a supported option; type ''? mesti'' for supported options."))
elseif lowercase(syst.zBC) == "bloch"
throw(ArgumentError("syst.zBC = \"Bloch\" but syst.kz_B is not given."))
end
zBC = syst.zBC
end
# Convert BC to take care of lowercase or uppercase
yBC = convert_BC(yBC, "y")
zBC = convert_BC(zBC, "z")
if ~use_2D_TM
xBC = convert_BC(xBC, "x")
# Check number of grid points with the boundary conditions
if nx_Ey != nx_Ez; throw(ArgumentError("Number of grids along x provided by syst.epsilon_yy and syst.epsilon_zz should be same.")); end
if ny_Ex != ny_Ez; throw(ArgumentError("Number of grids along y provided by syst.epsilon_xx and syst.epsilon_zz should be same.")); end
if nz_Ex != nz_Ey; throw(ArgumentError("Number of grids along z provided by syst.epsilon_xx and syst.epsilon_yy should be same.")); end
check_BC_and_grid(xBC, nx_Ex, nx_Ey, nx_Ez, "x")
check_BC_and_grid(yBC, ny_Ex, ny_Ey, ny_Ez, "y")
check_BC_and_grid(zBC, nz_Ex, nz_Ey, nz_Ez, "z")
if (isdefined(syst, :epsilon_xy) && ~isa(syst.epsilon_xy, Nothing) && ~(size(syst.epsilon_xy) == (nx_Ez, ny_Ez, nz_Ex)))
throw(ArgumentError("The size of syst.epsilon_xy should be should be (size(syst.epsilon_zz, 1), size(syst.epsilon_zz, 2), size(syst.epsilon_xx, 3)) = ($(size(syst.epsilon_zz, 1)), $(size(syst.epsilon_zz, 2)), $(size(syst.epsilon_xx, 3)))."))
end
if (isdefined(syst, :epsilon_xz) && ~isa(syst.epsilon_xz, Nothing) && ~(size(syst.epsilon_xz) == (nx_Ey, ny_Ex, nz_Ey)))
throw(ArgumentError("The size of syst.epsilon_xz should be should be (size(syst.epsilon_yy, 1), size(syst.epsilon_xx, 2), size(syst.epsilon_yy, 3)) = ($(size(syst.epsilon_yy, 1)), $(size(syst.epsilon_xx, 2)), $(size(syst.epsilon_yy, 3)))."))
end
if (isdefined(syst, :epsilon_yx) && ~isa(syst.epsilon_yx, Nothing) && ~(size(syst.epsilon_yx) == (nx_Ez, ny_Ez, nz_Ey)))
throw(ArgumentError("The size of syst.epsilon_yx should be should be (size(syst.epsilon_zz, 1), size(syst.epsilon_zz, 2), size(syst.epsilon_yy, 3)) = ($(size(syst.epsilon_zz, 1)), $(size(syst.epsilon_zz, 2)), $(size(syst.epsilon_yy, 3)))."))
end
if (isdefined(syst, :epsilon_yz) && ~isa(syst.epsilon_yz, Nothing) && ~(size(syst.epsilon_yz) == (nx_Ey, ny_Ex, nz_Ex)))
throw(ArgumentError("The size of syst.epsilon_yz should be should be (size(syst.epsilon_yy, 1), size(syst.epsilon_xx, 2), size(syst.epsilon_xx, 3)) = ($(size(syst.epsilon_yy, 1)), $(size(syst.epsilon_xx, 2)), $(size(syst.epsilon_xx, 3)))."))
end
if (isdefined(syst, :epsilon_zx) && ~isa(syst.epsilon_zx, Nothing) && ~(size(syst.epsilon_zx) == (nx_Ey, ny_Ez, nz_Ey)))
throw(ArgumentError("The size of syst.epsilon_zx should be should be (size(syst.epsilon_yy, 1), size(syst.epsilon_zz, 2), size(syst.epsilon_yy, 3)) = ($(size(syst.epsilon_yy, 1)), $(size(syst.epsilon_zz, 2)), $(size(syst.epsilon_yy, 3)))."))
end
if (isdefined(syst, :epsilon_zy) && ~isa(syst.epsilon_zy, Nothing) && ~(size(syst.epsilon_zy) == (nx_Ez, ny_Ex, nz_Ex)))
throw(ArgumentError("The size of syst.epsilon_zy should be should be (size(syst.epsilon_zz, 1), size(syst.epsilon_xx, 2), size(syst.epsilon_xx, 3)) = ($(size(syst.epsilon_zz, 1)), $(size(syst.epsilon_xx, 2)), $(size(syst.epsilon_xx, 3)))."))
end
end
# Defaults to no PML anywhere
if ~isdefined(syst, :PML)
pml = PML(0)
pml.direction = "all"
syst.PML = [pml]
elseif isa(syst.PML, PML)
# If user specifies PML structure instead of vector of PML strcture, transform it to vector of PML strcture.
syst.PML = [syst.PML]
end
# Parse the user-specified PML parameters to PML on the six sides
# PML_list = [xPML_low, xPML_high, yPML_low, yPML_high, zPML_low, zPML_high]
PML_list = Vector{PML}(undef, 6)
str_sides = ["-x", "+x", "-y", "+y", "-z", "+z"]
use_PML = false
for ii = 1:length(syst.PML)
PML_ii = syst.PML[ii]
# Check that the user did not accidentally use options only in mesti2s()
if isdefined(PML_ii, :npixels_spacer) && ~isa(PML_ii.npixels_spacer, Nothing)
@warn "syst.PML[$(ii)] field \"npixels_spacer\" is not used in mesti() and will be ignored."
PML_ii.npixels_spacer = nothing
end
# Number of PML pixels must be given
# Other PML parameters are optional and will be checked in mesti_build_fdfd_matrix()
if ~isdefined(PML_ii, :npixels)
throw(ArgumentError("syst.PML[$(ii)] must contain field \"npixels\"."))
end
if PML_ii.npixels != 0
use_PML = true
end
# Direction(s) where PML is placed must be given
if ~isdefined(PML_ii, :direction)
throw(ArgumentError("syst.PML[$(ii)] must contain field \"direction\"."))
elseif ~(lowercase(PML_ii.direction) in ["all", "x", "y", "z", "xy", "xz", "yz", "yx", "zx", "zy"])
throw(ArgumentError("syst.PML[$(ii)].direction = \"$(PML_ii.direction)\" is not a supported option; use \"all\", \"x\", \"y\", \"z\", \"xy\", \"xz\", \"yz\", \"yx\", \"zx\", or \"zy\"."))
end
# If PML is specified, we put it on both sides by default
if ~isdefined(PML_ii, :side)
PML_ii.side = "both"
elseif ~(lowercase(PML_ii.side) in ["both", "-", "+"])
throw(ArgumentError("syst.PML[$(ii)].side = \"$(PML_ii.side)\" is not a supported option; use \"both\", \"-\", or \"+\"."))
end
# Convert [PML_ii.direction and PML_ii.side] to a list of the PML locations
# 1=xPML_low, 2=xPML_high, 3=yPML_low, 4=yPML_high, 5=zPML_low, 6=zPML_high
if lowercase(PML_ii.direction) == "all" # x & y & z
if lowercase(PML_ii.side) == "both" # low & high
if ~use_2D_TM
ind_ii = [1,2,3,4,5,6]
else
# For 2D TM fields Ex(y,z), only put PML on y and z directions
ind_ii = [3,4,5,6]
end
elseif PML_ii.side == "-"
if ~use_2D_TM
ind_ii = [1,3,5]
else
# For 2D TM fields Ex(y,z), only put PML on y and z directions
ind_ii = [3,5]
end
else # PML_ii.side == "+"
if ~use_2D_TM
ind_ii = [2,4,6]
else
# For 2D TM fields Ex(y,z), only put PML on y and z directions
ind_ii = [4,6]
end
end
elseif lowercase(PML_ii.direction) == "x"
if use_2D_TM
@warn "syst.PML[$(ii)].direction = \"x\" would not apply to 2D cases; will be ignored."
end
if lowercase(PML_ii.side) == "both" # low & high
ind_ii = [1,2]
elseif PML_ii.side == "-"
ind_ii = 1
else # PML_ii.side = "+"
ind_ii = 2
end
elseif lowercase(PML_ii.direction) == "y"
if lowercase(PML_ii.side) == "both" # low & high
ind_ii = [3,4]
elseif lowercase(PML_ii.side) == "-"
ind_ii = 3
else # PML_ii.side = "+"
ind_ii = 4
end
elseif lowercase(PML_ii.direction) == "z"
if lowercase(PML_ii.side) == "both" # low & high
ind_ii = [5,6]
elseif lowercase(PML_ii.side) == "-"
ind_ii = 5
else # PML_ii.side = "+"
ind_ii = 6
end
elseif lowercase(PML_ii.direction) == "xy" || lowercase(PML_ii.direction) == "yx"
if use_2D_TM
@warn "syst.PML[$(ii)].direction = \"$(PML_ii.side)\"; in 2D cases, PML in x-direction will not be applied and ignored."
end
if lowercase(PML_ii.side) == "both" # low & high
ind_ii = [1,2,3,4]
elseif lowercase(PML_ii.side) == "-"
ind_ii = [1,3]
else # PML_ii.side = "+"
ind_ii = [2,4]
end
elseif lowercase(PML_ii.direction) == "xz" || lowercase(PML_ii.direction) == "zx"
if use_2D_TM
@warn "syst.PML[$(ii)].direction = \"$(PML_ii.side)\"; in 2D cases, PML in x-direction will not be applied and ignored."
end
if lowercase(PML_ii.side) == "both" # low & high
ind_ii = [1,2,5,6]
elseif lowercase(PML_ii.side) == "-"
ind_ii = [1,5]
else # PML_ii.side = "+"
ind_ii = [2,6]
end
else # lowercase(PML_ii.direction) == "yz" || lowercase(PML_ii.direction) == "zy"
if lowercase(PML_ii.side) == "both" # low & high
ind_ii = [3,4,5,6]
elseif lowercase(PML_ii.side) == "-"
ind_ii = [3,5]
else # PML_ii.side = "+"
ind_ii = [4,6]
end
end
# Specify PML at those locations
for jj = 1:length(ind_ii)
ind_side = ind_ii[jj]
# Check that PML has not been specified at that location yet
if isassigned(PML_list,ind_side)
throw(ArgumentError("PML on $(str_sides[ind_side]) side is specified more than once in syst.PML."))
end
PML_list[ind_side] = deepcopy(PML_ii)
end
end
for ii = 1:length(PML_list)
if ~isassigned(PML_list, ii)
PML_list[ii] = PML(0)
end
PML_list[ii].side = string(str_sides[ii][1])