-
Notifications
You must be signed in to change notification settings - Fork 3
/
mesti2s.jl
2803 lines (2616 loc) · 174 KB
/
mesti2s.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 channel_type
export channel_index
export wavefront
# Export a function mesti2s()
export mesti2s
mutable struct channel_type
# A composite data type to specify channel type
# See also: mesti and mesti2s
side::String
polarization::String
channel_type() = new()
end
mutable struct channel_index
# A composite data type to specify channel index
# See also: mesti and mesti2s
# Used in 3D systems
ind_low_s::Vector{Int64}
ind_low_p::Vector{Int64}
ind_high_s::Vector{Int64}
ind_high_p::Vector{Int64}
# Used in 2D systems
ind_low::Vector{Int64}
ind_high::Vector{Int64}
channel_index() = new()
end
mutable struct wavefront
# A composite data type to specify wavefront, which is linear combinations of channels
# See also: mesti and mesti2s
# Used in 3D systems
v_low_s::Union{Array{Int64,2}, Array{Float64,2}, Array{ComplexF64,2}}
v_low_p::Union{Array{Int64,2}, Array{Float64,2}, Array{ComplexF64,2}}
v_high_s::Union{Array{Int64,2}, Array{Float64,2}, Array{ComplexF64,2}}
v_high_p::Union{Array{Int64,2}, Array{Float64,2}, Array{ComplexF64,2}}
# Used in 2D systems
v_low::Union{Array{Int64,2}, Array{Float64,2}, Array{ComplexF64,2}}
v_high::Union{Array{Int64,2}, Array{Float64,2}, Array{ComplexF64,2}}
wavefront() = new()
end
"""
MESTI2S Solves frequency-domain scattering problems in a two-sided geometry.
---3D field profile---
(Ex, Ey, Ez, channels, info) = MESTI2S(syst, input) 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". epsilon(x,y,z) must have homogeneous spaces on the low (-z)
and high (+z) sides, with an outgoing boundary in z for the scattered waves and a closed
boundary in x and y.
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.
The incident wavefronts from low and/or high are specified by variable "input".
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".
---2D TM field profile---
(Ex, channels, info) = MESTI2S(syst, input) returns the spatial
field profiles of Ex(y,z) for scattering problems of 2D transverse-magnetic (TM) fields satisfying
[- (d/dy)^2 - (d/dz)^2 - (omega/c)^2*epsilon_xx(y,z)] Ex(y,z) = 0,
where Ex is each a superposition of incident and scattered waves.
The returned "Ex" is a 3D array, with Ex(:,:,i) being the total (incident + scattered) field profile
of Ex given the i-th incident wavefront.
---Scattering matrix S---
(S, channels, info) = MESTI2S(syst, input, output) returns the scattering matrix
S, where "input" and "output" specify either the list of input/output channels or
the input/output wavefronts. When the MUMPS solver is available,
this is typically done by computing the Schur complement of an augmented
matrix K through a partial factorization.
(Ex, Ey, Ez, channels, info) = MESTI2S(syst, input, opts),
(Ex, channels, info) = MESTI2S(syst, input, opts), and
(S, channels, info) = MESTI2S(syst, input, output, opts) allow detailed options to
be specified with structure "opts".
In mesti2s(), for 3D cases, the boundary condition in x and y must be closed
e.g., periodic or PEC. For 2D cases, the boundary condition in y must be closed.
Given the closed boundary, the set of transverse modes forms a
complete and orthonormal basis of propagating and evanescent channels.
The inputs and outputs are specified in the basis of these propagating
channels, with coefficients normalized with respect to the flux in the
longitudinal (z) direction. Properties of those channels are given by
mesti_build_channels().
When in 3D cases an open boundary in x or y is of interest (in 2D cases an open boundary
in y is of interest), the appropriate input/output basis depends on the problem,
so the user should use the more general function mesti() and will need to specify
the input and output matrices B and C.
MESTI only considers nonmagnetic materials.
This file builds the input and output channels using mesti_build_channels(),
builds the matrices B and C, and then calls mesti() to solve the scattering
problems.
=== 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):
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.
Outside the scattering region (with z < 0 and z > H), epsilon_xx(x,y,z)
is given by scalars syst.epsilon_low and syst.epsilon_high for the
homogeneous spaces on the two sides. Note that nx_Ez = 0 corresponds
to H = 0 (ie, no scattering region) and is allowed.
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.epsilon_low (real scalar; required):
Relative permittivity of the homogeneous space on the low.
syst.epsilon_high (real scalar or nothing; optional):
Relative permittivity of the homogeneous space on the high. If
syst.epsilon_high is not given or is nothing, for 3D cases the system
will be one-sided, terminated on the high with a PEC boundary with
Ez(n,m,l) = 0 at l = nz_Ez + 1, and for 2D TM caes the system
will be one-sided, terminated on the high with a PEC boundary with
Ex(m,l) = 0 at l = nx_Ex + 1.
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, 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.xBC (string or nothing; required unless syst.kx_B is specified):
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.
Note that this xBC also defines a complete and orthonormal set of
transverse modes, upon which the input and output channels in
arguments "input" and "output" are defined; mesti2s() does not support PML
in x direction because a closed boundary is necessary for defining
such a transverse basis.
Here, syst.xBC is required, with no default choice (except when
syst.kx_B is given, in which case syst.xBC = "Bloch" is automatically
used).
For 2D TM fields, xBC = nothing.
syst.yBC (string; optional):
Boundary condition in y direction, analogous to syst.xBC.
syst.kx_B (real scalar; 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.
syst.ky_B (real scalar; optional):
Bloch wave number in y direction, analogous to syst.kx_B.
syst.zPML (a vector of PML structure; required):
Parameters of the perfectly matched layer (PML) used to simulate an
open boundary, which attenuates outgoing waves with minimal reflection.
Note that in mesti2s(), the PML is placed in the homogeneous spaces
specified by syst.epsilon_low and syst.epsilon_high, outside of the
scattering region specified by syst.epsilon_xx/Ey/Ez. (This is different
from the more general function mesti(), where syst.epsilon_xx/Ey/Ez
specifies the entire simulation domain, so PML is placed inside
syst.epsilon_xx/Ey/Ez.)
To use the same PML on both sides or to use PML on the low of a
one-sided geometry, set syst.zPML to be a scalar structure with the
following fields:
npixels (positive integer scalar; required): Number of PML pixels.
This number of pixels is added in addition to the
scattering region.
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.
When there is just one PML element in the vector,
syst.zPML = [PML],
then the same PML parameters are used in the both sides of two-sided geometry,
or the PML parameters are used in low side of one-sided geometry.
When multiple sets of PML parameters are used in the two-sided geometry
(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.zPML = [PML_low, PML_high],
with PML_low and PML_high each being a structure containing the above
fields; they can specify different PML parameters on low and high sides.
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.
input (channel_type, channel_index, or wavefront structure; required):
The set of input channels or input wavefronts
To specify all propagating channels on one side or on both sides for
polarizations, use ''input = channel_type()'', then it contains the following fields:
input.side (string, required): specify all propagating channels on
sides. Available choices are:
"low" - specify all (input.polarization)-polarization
propagating channels on the low side
"high" - specify all (input.polarization)-polarization
propagating channels on the high side
"both" - specify all (input.polarization)-polarization
propagating channels on the both sides
input.polarization (string, optional): specify polarizations on
sides for 3D systems. Available choices are:
"s" - specify s-polarization channels on the input.side
"p" - specify p-polarization channels on the input.side
"both" - specify both polarizations (s and p) channels on
on the input.side
By default, input.polarization = "both";
To specify a subset of the propagating channels use ''input = channel_index()'',
then it contains the following fields:
For 3D systems:
input.ind_low_s (integer vector): Vector containing the indices of
propagating channels incident on the low side with s-polarization
input.ind_low_p (integer vector): Vector containing the indices of
propagating channels incident on the low side with p-polarization
input.ind_high_s (integer vector): Vector containing the indices of
propagating channels incident on the high side with s-polarization
input.ind_high_p (integer vector): Vector containing the indices of
propagating channels incident on the high side with p-polarization
One can provide only one or more of them.
The user can first use mesti_build_channels() to get the indices, wave
numbers, and transverse profiles of the propagating channels; based on
that, the user can specify the list of channels of interest through
input.ind_low_s, input.ind_low_p, input.ind_high_s, and input.ind_high_p
above, or a list of customized wavefronts through input.v_low_s,
input.v_low_p, input.v_high_s, or input.v_high_p below.
For 2D TM case:
input.ind_low (integer vector): Vector containing the indices of
propagating channels incident on the low side.
input.ind_high (integer vector): Vector containing the indices of
propagating channels incident on the high side.
To specify a custom input wavefronts, a superposition of multiple
propagating channels, use ''input = wavefront()'', then it contains the following fields:
For 3D systems:
input.v_low_s (numeric matrix): Matrix where each column specifies the
coefficients of s-polarization propagating channels on the low
for one input wavefront from the low; the wavefront is a superposition
of all propagating channels with the superposition coefficients
given by that column of input.v_low_s. size(input.v_low_s, 1) must equal
N_prop_low, the total number of propagating channels on the low;
size(input.v_low_s, 2) is the number of input wavefronts.
input.v_low_p (numeric matrix): Analogous to to input.v_low_s, but
specifying input p-polarization wavefronts from the low instead.
input.v_high_s (numeric matrix): Analogous to to input.v_low_s, but
specifying input s-polarization wavefronts from the high instead.
input.v_high_p (numeric matrix): Analogous to to input.v_low_s, but
specifying input p-polarization wavefronts from the high instead.
Note that the input wavefronts from the low and the input wavefronts
from the high and two polarizations are treated as separate inputs. In other
words, each input either comes from the low or comes from the high on one
polarization. If an input with incidence from both sides or both polarization
is of interest, the user can manually superimpose results from the separate
computations.
For 2D TM case:
input.v_low (numeric matrix): Matrix where each column specifies the
coefficients of propagating channels on the low for one input
wavefront from the low, with the superposition coefficients given by that
column of input.v_low. size(input.v_low, 1) must equal N_prop_low, the total
number of propagating channels on the low; size(input.v_L, 2) is the number
of input wavefronts.
input.v_high (numeric matrix): Analogous to to input.v_low, but specifying
input wavefronts from the high instead.
output (channel_type, channel_index, wavefront structure, or nothing; optional):
The set of output channels or output wavefronts.
When out = nothing or when out is omitted, no output projection is used, and
the spatial field profiles Ex(x,y,z), Ey(x,y,z), and Ez(x,y,z) corresponding to
the set of inputs are returned.
When out is given, the scattering matrix is returned, with the output
basis of the scattering matrix specified by "output". In this case, out
follows the same syntax as the argument "input".
To specify all propagating channels on one side or on both sides for
polarizations, use ''output = channel_type()'', then it contains the following fields:
output.side (string, required): specify all propagating channels on sides.
Available choices are:
"low" - specify all propagating channels on the low side for
polarization specified by output.polarization
"high" - specify all propagating channels on the high side for
polarization specified by output.polarization
"both" - specify all propagating channels on the both sides
(low and high) for polarization specified by
output.polarization
output.polarization (string, optional): specify polarizations on sides for 3D systems.
Available choices are:
"s" - specify s-polarization channels on the output.side
"p" - specify p-polarization channels on the output.side
"both" - specify both polarizations (s and p) channels on
on the output.side
By default, output.polarization = "both";
To specify a subset of the propagating channels use ''output = channel_index()'',
then it contains the following fields:
For 3D systems:
output.ind_low_s (integer vector): Vector containing the indices of
propagating channels incident on the low side with s-polarization
output.ind_low_p (integer vector): Vector containing the indices of
propagating channels incident on the low side with p-polarization
output.ind_high_s (integer vector): Vector containing the indices of
propagating channels incident on the high side with s-polarization
output.ind_high_p (integer vector): Vector containing the indices of
propagating channels incident on the high side with p-polarization
One can provide only one or more of them.
The user can first use mesti_build_channels() to get the indices, wave
numbers, and transverse profiles of the propagating channels; based on
that, the user can specify the list of channels of interest through
output.ind_low_s, output.ind_low_p, output.ind_high_s, and output.ind_high_p
above, or a list of customized wavefronts through output.v_low_s, output.v_low_p,
output.v_high_s, or output.v_high_p below.
For 2D TM case:
output.ind_low (integer vector): Vector containing the indices of
propagating channels incident on the low side.
output.ind_high (integer vector): Vector containing the indices of
propagating channels incident on the high side.
To specify a custom output wavefronts, a superposition of multiple
propagating channels, use ''output = wavefront()'', then it contains the following fields:
For 3D systems:
output.v_low_s (numeric matrix): Matrix where each column specifies the
coefficients of s-polarization propagating channels on the low
for one output wavefront from the low; the wavefront is a superposition
of all propagating channels with the superposition coefficients
given by that column of output.v_low_s. size(output.v_low_s, 1) must equal
N_prop_low, the total number of propagating channels on the low;
size(output.v_low_s, 2) is the number of output wavefronts.
output.v_low_p (numeric matrix): Analogous to to output.v_low_s, but
specifying output p-polarization wavefronts from the low instead.
output.v_high_s (numeric matrix): Analogous to to output.v_low_s, but
specifying output s-polarization wavefronts from the high instead.
output.v_high_p (numeric matrix): Analogous to to output.v_low_s, but
specifying output p-polarization wavefronts from the high instead.
Here, each row of the scattering matrix corresponds to projection onto an
output wavefront specified by a column of output.v_low_s, output.v_low_p,
output.v_high_s, or output.v_high_p.
For 2D TM case:
input.v_low (numeric matrix): Matrix where each column specifies the
coefficients of propagating channels on the low for one input
wavefront from the low, with the superposition coefficients given by that column of input.v_low.
size(input.v_low, 1) must equal N_prop_low, the total number of
propagating channels on the low; size(input.v_L, 2) is the number
of input wavefronts.
input.v_high (numeric matrix): Analogous to to input.v_low, but specifying
input wavefronts from the high instead.
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.nz_low (non-negative integer scalar; optional, defaults to 0):
Number of pixels of homogeneous space on the low (syst.epsilon_low) to
include when returning the spatial field profile; not used for
scattering matrix computations.
opts.nz_high (non-negative integer scalar; optional, defaults to 0):
Number of pixels of homogeneous space on the high (syst.epsilon_high) to
include when returning the spatial field profile; not used for
scattering matrix computations. Note that opts.nx_high can still be used
in one-sided geometries where syst.epsilon_high is not given; the field
profile on the high is simply zero in such case.
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.symmetrize_K (boolean scalar; optional):
Whether or not to pad input and/or output channels and perform
permutations to make matrix K = [A,B;C,0] symmetric when computing its
Schur complement, which lowers computing time and memory usage. Such
channel padding and permutation is reversed afterwards and does not
affect what mesti2s() returns.
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().
opts.clear_memory (boolean scalar; optional, defaults to true):
Whether or not to clear variables inside mesti() to reduce peak memory
usage.
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_continuous_dispersion (boolean scalar; optional):
Whether to use the dispersion equation of the continuous wave equation
when building the input/output channels. Defaults to false, in which
case the finite-difference dispersion is used.
opts.n0 (real numeric scalar; optional, defaults to 0):
Center of the 1D transverse mode profile with periodic or Bloch periodic
boundary condition, u(m,a) = exp(i*kx(a)*syst.dx*(m-m0))/sqrt(nx),
where kx(a) = syst.kx_B + a*(2*pi/nx*syst.dx) and nx = nx_Ex = nx_Ey = nx_Ez.
opts.m0 (real numeric scalar; optional, defaults to 0):
Center of the 1D transverse mode profile with periodic or Bloch periodic
boundary conditions, analogous to opts.n0.
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 out---
field_profiles (4D array):
For field-profile computations (i.e., when "out" is not given), the
returned field_profiles is a 4D array containing the field profiles, such
that field_profiles(:,:,:,a) is the total (incident + scattered) field of
Ex/Ey/Ez corresponding to the a-th input wavefront, with
Recall that nz_Ez = nz_Ex + 1 = nz_Ey + 1 = H/syst.dx + 1 for a
two-sided geometry; nz_Ez = nz_Ex = nz_Ey = H/syst.dx for a one-sided geometry.
By default, opts.nz_low = opts.nz_high = 0, and field_profiles only contain
field_profiles in the scattering region. When opts.nz_low and opts.nz_high are
specified, field_profiles will also contain opts.nz_low pixels in the homogeneous
space on the low, opts.nz_high pixel in the homogeneous space on the high.
To be specific, the field_profiles contains:
For 3D systems:
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 Ex component
(nx_Ey, ny_Ey, nz_Ey, number of inputs) = size(Ey)
Ez (4D array):
electrical field profile for Ex 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)
---or when users do not specify out---
S (numeric matrix):
For scattering-matrix computations (i.e., when ''output'' is given), S is the
scattering matrix, such that S(b,a) is the flux-normalized field
coefficient in the b-th propagating output channel (or the b-th output
wavefront) given incident wave in the a-th propagating input channel (or
the a-th input wavefront).
When all propagating channels on one side or on both sides are
requested, e.g. with ''input = channel_type()'', ''input.side = high'' or
''output = channel_type()'', ''output.side = both'', matrix S includes
channels.low.N_prop propagating channels on the low and/or
2*channels.high.N_prop on the high, with longitudinal wave
numbers given by vectors channels.low.kzdx_prop and channels.high.kzdx_prop.
The phases of the elements of the scattering matrix depend on the
reference plane. For channels on the low, the reference plane is at
z = 0. For channels on the high, the reference plane is at z = L.
When a subset of the propagating channels are requested, in 3D
''input = channel_index()'' and ''output = channel_index()''
with in.ind_low_s, in.ind_low_p, in.ind_high_s, in.ind_high_p,
out.ind_low_s, out.ind_low_p, out.ind_low_p and/or out.ind_high_p,
matrix S includes such subset of the propagating channels.
In 2D, it used with in.ind_low, in.ind_high, out.ind_low, and out.ind_high.
When the input wavefronts and/or output basis, i.e. ''input = wavefront()''
and/or ''output = wavefront()'' is specified in 3D by input.v_low_s, input.v_low_p,
input.v_high_s, input.v_high_p, output.v_low_s, output.v_low_p, output.v_high_s,
and/or output.v_high_p. In 2D, it used with input.v_low, input.v_high, output.v_low,
and output.v_high. Matrix S is the full scattering matrix with the input and/or output
basis changed.
channels (Channel structure):
A structure returned by function mesti_build_channels() that contains
properties of the propagating and evanescent channels of the homogeneous
spaces on the low and high. Type "? mesti_build_channels" for more
information.
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):
info.timing_build (non-negative scalar):
info.timing_analyze (non-negative scalar):
info.timing_factorize (non-negative scalar):
info.timing_total (non-negative scalar):
info.zPML (two-element cell array);
PML parameters on the sides of z direction.
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_channels, mesti2s
"""
function mesti2s(syst::Syst, input::Union{channel_type, channel_index, wavefront}, output::Union{channel_type, channel_index, wavefront, Nothing}, opts::Union{Opts, Nothing})
# Make deepcopy of them to avoid mutating input argument
syst = deepcopy(syst); input = deepcopy(input); output = deepcopy(output); opts = deepcopy(opts)
## Part 1.1: Check validity of syst, assign default values to its fields, and parse BC and PML specifications
t0 = time()
if ~(isdefined(syst, :epsilon_low) && isa(syst.epsilon_low, Real)); throw(ArgumentError("Input argument syst must have field \"epsilon_low\" and be a real scalar.")); end
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 && ~(isdefined(syst, :epsilon_xx) && isdefined(syst, :epsilon_yy) && isdefined(syst, :epsilon_zz))
throw(ArgumentError("Input argument syst must have field \"epsilon_xx\", \"epsilon_yy\", and \"epsilon_zz\" for 3D systems."))
end
# Check that the user did not accidentally use options only in mesti()
if isdefined(syst, :PML)
throw(ArgumentError("syst.PML is not used in mesti2s(); use syst.zPML instead."))
elseif isdefined(syst, :zBC)
throw(ArgumentError("syst.zBC is not supported in mesti2s(); use mesti() if the boundary condition behind PML is not PEC."))
elseif isdefined(syst, :kz_B)
throw(ArgumentError("syst.kz_B is not supported in mesti2s(); use mesti() if Bloch periodic BC in z is needed."))
elseif isdefined(syst, :PML_type) && (lowercase(syst.PML_type) != "upml")
throw(ArgumentError("syst.PML_type = \"UPML\" is the only optiona in mesti2s(); use mesti() if SC-PML is needed."))
end
# syst.epsilon_high is an optional argument; determines whether the system is one-sided
if ~isdefined(syst, :epsilon_high) || isa(syst.epsilon_high, Nothing)
two_sided = false
syst.epsilon_high = nothing
elseif ~isa(syst.epsilon_high, Real)
throw(ArgumentError("syst.epsilon_high must be a real scalar, if given."))
else
two_sided = true
end
if ~use_2D_TM
# 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)
if maximum([nx_Ex*ny_Ex, nx_Ey*ny_Ey, nx_Ez*ny_Ez]) == 0
throw(ArgumentError("Total number of pixel cannot be 0 for all Ex, Ey, and Ez components in the tranverse xy plane."))
end
else
# Number of grid points in y and z for Ex
(ny_Ex, nz_Ex) = size(syst.epsilon_xx)
end
dn = 0.5 # the source/detection is half a pixel away from z = 0 and z = L
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_channels() uses kx_B*periodicity as the input arguments xBC for Bloch BC
xBC = (syst.kx_B)*(nx_Ex*syst.dx) # dimensionless
else
if ~isdefined(syst, :xBC)
throw(ArgumentError("Input argument syst must have non-empty field \"xBC\" when syst.kx_B is not given."))
elseif ~(lowercase(syst.xBC) in ["bloch", "periodic", "pec", "pmc", "pecpmc", "pmcpec"])
throw(ArgumentError("syst.xBC = \"$(syst.xBC)\" is not a supported option; type ''? mesti2s'' 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_channels() uses ky_B*periodicity as the input arguments yBC 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
if ~isdefined(syst, :yBC)
throw(ArgumentError("Input argument syst must have non-empty field \"yBC\" when syst.ky_B is not given."))
elseif ~(lowercase(syst.yBC) in ["bloch", "periodic", "pec", "pmc", "pecpmc", "pmcpec"])
throw(ArgumentError("syst.yBC = \"$(syst.yBC)\" is not a supported option; type ''? mesti2s'' 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 zPML
if ~isdefined(syst, :zPML)
throw(ArgumentError("syst.zPML must be given."))
elseif isa(syst.zPML, Nothing)
throw(ArgumentError("syst.zPML must be a PML structure or a vector of PML structure."))
elseif isa(syst.zPML, PML)
syst.zPML = [syst.zPML]
end
if length(syst.zPML) > 2 || length(syst.zPML) == 0
throw(ArgumentError("The length of syst.zPML must be 1 or 2."))
end
if ~(length(syst.zPML)==2) && two_sided
if isdefined(syst.zPML, :side) && syst.zPML.side != "both"
throw(ArgumentError("syst.zPML.side should be assigned as \"both\", if given, when the user just provided only one PML object is provided for a two-sided geometry."))
end
# Apply the same zPML on both sides
syst.zPML = [deepcopy(syst.zPML[1]), deepcopy(syst.zPML[1])]
syst.zPML[1].side = "-"; syst.zPML[2].side = "+"
elseif length(syst.zPML)==2 && ~two_sided
throw(ArgumentError("For a one-sided geometry, the boundary condition on the high side is PEC; syst.zPML only specifies PML on the low side and must be a PML structure or a vector containing only one PML structure."))
end
# Convert BC to take care of lowercase or uppercase
if ~use_2D_TM
xBC = convert_BC(xBC, "x")
end
yBC = convert_BC(yBC, "y")
if ~use_2D_TM
# 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("PEC", 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
# nz_extra is the number of homogeneous-space pixels to be added to syst.epsilon_xx, epsilon_yx, epsilon_zx, epsilon_xy, syst.epsilon_yy, epsilon_zy, epsilon_xz, epsilon_yz and syst.epsilon_zz in z direction.
# The two elements of nz_extra corresponds to the low and high sides.
if two_sided
n_sides = 2
nz_extra = Vector{Int64}(undef, 2)
str_zBC = ["PML", "PML"] # to be used for printing later
else
n_sides = 1
nz_extra = Vector{Int64}(undef, 2)
nz_extra[2] = 0
str_zBC = ["PML", "PEC"]
end
syst.PML = Vector{PML}() # to be used in mesti()
# Loop over syst.zPML and handle PML parameters, if specified
str_sides = ["-", "+"]
for ii = 1:n_sides
PML_ii = syst.zPML[ii]
if isdefined(PML_ii, :direction) && PML_ii.direction != "z"
throw(ArgumentError("syst.zPML[$(ii)].direction must be \"z\", if given."))
else
PML_ii.direction = "z" # to be used in mesti()
end
if isdefined(PML_ii, :side) && PML_ii.side != str_sides[ii]
throw(ArgumentError("syst.zPML[$(ii)].side must be $(str_sides[ii]), if given."))
else
PML_ii.side = str_sides[ii] # to be used in mesti()
end
# Number of PML pixels must be given
# Other fields are optional and will be checked in mesti_build_fdfd_matrix()
if ~isdefined(PML_ii, :npixels)
throw(ArgumentError("syst.zPML[$(ii)] must contain field \"npixels\"."))
end
if ~(PML_ii.npixels >= 0)
throw(ArgumentError("syst.zPML[$(ii)].npixels must be a non-negative integer scalar."))
end
# Defaults to no spacer
if ~isdefined(PML_ii, :npixels_spacer)
PML_ii.npixels_spacer = 0
end
if ~(PML_ii.npixels_spacer >= 0)
throw(ArgumentError("syst.zPML[$(ii)].npixels_spacer must be a non-negative integer scalar, if given."))
end
# number of pixels in z to be added outside of syst.epsilon_xx, syst.epsilon_yy, and syst.epsilon_zz
nz_extra[ii] = 1 + PML_ii.npixels + PML_ii.npixels_spacer
PML_ii.npixels_spacer = nothing # these will not be used in mesti()
push!(syst.PML, PML_ii) # add this PML structure to the vector of PML structure
end
syst.zPML = nothing # mesti() will throw warning if syst.zPML is given
if ~use_2D_TM
# Total number of pixels in z
nz_tot_Ex = nz_Ex + sum(nz_extra)
nz_tot_Ey = nz_Ey + sum(nz_extra)
nz_tot_Ez = nz_Ez + sum(nz_extra)
# Total number of grid points for Ex, Ey, and Ez
nt_tot_Ex = nx_Ex*ny_Ex*nz_tot_Ex
nt_tot_Ey = nx_Ey*ny_Ey*nz_tot_Ey
nt_tot_Ez = nx_Ez*ny_Ez*nz_tot_Ez
else
# Total number of pixels in z
nz_tot_Ex = nz_Ex + sum(nz_extra)
# Total number of grid points for Ex
nt_tot_Ex = ny_Ex*nz_tot_Ex
end
## Part 1.2: Check the arguments opts and assign default values
## We continue from here next time
if opts == nothing
opts = Opts()
end
# Check that the user did not accidentally use options only in mesti()
if isdefined(opts, :prefactor) && ~isa(opts.prefactor, Nothing)
throw(ArgumentError("opts.prefactor is not used in mesti2s(); the -2i prefactor is automatically included."))
end
# We return the scattering matrix S=C*inv(A)*B if output is given from input argument; else we return the spatial field profiles
# This opts.return_field_profile is not specified by the user; it will be returned as info.opts.return_field_profile to help debugging
opts.return_field_profile = isa(output, Nothing)
# By default, for field-profile computation, we only return result within the scattering region, setting nz_low = nz_high = 0.
if opts.return_field_profile
if ~isdefined(opts, :nz_low) || isa(opts.nz_low, Nothing)
opts.nz_low = 0
elseif ~(opts.nz_low >= 0)
throw(ArgumentError("opts.nz_low must be a non-negative integer scalar, if given."))
end
if ~isdefined(opts, :nz_high) || isa(opts.nz_high, Nothing)
opts.nz_high = 0
elseif ~(opts.nz_high >= 0 )
throw(ArgumentError("opts.nz_high must be a non-negative integer scalar, if given."))
end
else
if isdefined(opts, :nz_low) && ~isa(opts.nz_low, Nothing)
@warn "opts.nz_low is not used for scattering matrix computation; will be ignored."
opts.nz_low = nothing
end
if isdefined(opts, :nz_high) && ~isa(opts.nz_high, Nothing)
@warn ("opts.nz_high is not used for scattering matrix computation; will be ignored.")
opts.nz_high = nothing
end
end
# Turn on verbal output by default
if ~isdefined(opts, :verbal)
opts.verbal = true
elseif ~isa(opts.verbal, Bool)
throw(ArgumentError("opts.verbal must be a boolean, if given."))
end
# By default, we don't clear syst in the caller's workspace
if ~isdefined(opts, :clear_syst)
opts.clear_syst = false
elseif ~isa(opts.clear_syst, Bool)
throw(ArgumentError("opts.clear_syst must be a boolean, if given."))
end
# By default, we will clear internal variables
if ~isdefined(opts, :clear_memory)
opts.clear_memory = true
elseif ~isa(opts.clear_memory, Bool)
throw(ArgumentError("opts.clear_memory must be a boolean, if given."))
end
# Use MUMPS for opts.solver when it is available
MUMPS_available = haskey(ENV,"MUMPS_PREFIX")
solver_specified = true
if ~isdefined(opts, :solver) || isa(opts.solver, Nothing)
solver_specified = false
if MUMPS_available
opts.solver = "MUMPS"
else
opts.solver = "JULIA"
end
else
opts.solver = uppercase(opts.solver)
if ~(opts.solver in ["MUMPS", "JULIA"])
throw(ArgumentError("opts.solver = \"$(opts.solver)\" is not a supported option; use \"MUMPS\" or \"JULIA\"."))
elseif opts.solver == "MUMPS" && ~MUMPS_available
throw(ArgumentError("opts.solver = \"$(opts.solver)\" but the Julia environment variable \"MUMPS_PREFIX\" is not specified."))
end
end
method_specified = true
if ~isdefined(opts, :method) || isa(opts.method, Nothing)
method_specified = false
# By default, if the argument ''output'' is not given or if
# opts.iterative_refinement = true, then "factorize_and_solve" is used.
# Otherwise, then "APF" is used when opts.solver = "MUMPS", and
# "C*inv(U)*inv(L)*B" is used when opts.solver = "JULIA".
if opts.return_field_profile || (isdefined(opts, :iterative_refinement) && opts.iterative_refinement == true)
opts.method = "factorize_and_solve"
elseif isdefined(opts, :store_ordering) && opts.store_ordering == true
opts.method = "APF"
elseif isdefined(opts, :analysis_only) && opts.analysis_only == true
opts.method = "APF"
else
if opts.solver == "MUMPS"