-
Notifications
You must be signed in to change notification settings - Fork 587
/
meep.hpp
2364 lines (2027 loc) · 99 KB
/
meep.hpp
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
/* Copyright (C) 2005-2021 Massachusetts Institute of Technology
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2, or (at your option)
% any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software Foundation,
% Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/
#ifndef MEEP_H
#define MEEP_H
#include <functional>
#include <limits>
#include <memory>
#include <unordered_map>
#include <vector>
#include <stdio.h>
#include <stddef.h>
#include <math.h>
#include "meep/vec.hpp"
#include "meep/mympi.hpp"
#include "meep/meep-config.h"
namespace meep {
/* The (time-domain) fields arrays of the fields_chunk class as well
as the material arrays of the structure_chunk class chi1inv,
chi3, sigma, etc. can be stored using single-precision floating
point rather than double precision (the default). The reduced
precision can provide for up to a factor of 2X improvement in the
time-stepping rate with generally negligible loss in accuracy. */
#if MEEP_SINGLE // set to 1 via configure --enable-single
typedef float realnum;
#else
typedef double realnum;
#endif
#define MEEP_MIN_OUTPUT_TIME 4.0 // output no more often than this many seconds
extern int
verbosity; // if 0, suppress all non-error messages from Meep; 1 is default, 2 is debug output
const double pi = 3.141592653589793238462643383276;
const double infinity = HUGE_VAL;
#ifdef NAN
const double nan = NAN;
#else
const double nan = -7.0415659787563146e103; // ideally, a value never encountered in practice
#endif
class h5file;
// Defined in monitor.cpp
void matrix_invert(std::complex<double> (&Vinv)[9], std::complex<double> (&V)[9]);
// Defined in dft.cpp
std::vector<double> linspace(double freq_min, double freq_max, size_t Nfreq);
double pml_quadratic_profile(double, void *);
/* generic base class, only used by subclassing: represents susceptibility
polarizability vector P = chi(omega) W (where W = E or H). */
class susceptibility {
public:
susceptibility() {
id = cur_id++;
ntot = 0;
next = NULL;
FOR_COMPONENTS(c) FOR_DIRECTIONS(d) {
sigma[c][d] = NULL;
trivial_sigma[c][d] = true;
}
}
susceptibility(const susceptibility &s) {
id = s.id;
ntot = s.ntot;
next = NULL;
FOR_COMPONENTS(c) FOR_DIRECTIONS(d) {
sigma[c][d] = NULL;
trivial_sigma[c][d] = true;
}
}
virtual susceptibility *clone() const;
virtual ~susceptibility() {
FOR_COMPONENTS(c) FOR_DIRECTIONS(d) { delete[] sigma[c][d]; }
delete next;
}
int get_id() const { return id; }
bool operator==(const susceptibility &s) const { return id == s.id; };
// Returns the 1st order nonlinear susceptibility (generic)
virtual std::complex<realnum> chi1(realnum freq, realnum sigma = 1);
// update all of the internal polarization state given the W field
// at the current time step, possibly the previous field W_prev, etc.
virtual void update_P(realnum *W[NUM_FIELD_COMPONENTS][2],
realnum *W_prev[NUM_FIELD_COMPONENTS][2], realnum dt, const grid_volume &gv,
void *P_internal_data) const {
(void)P;
(void)W;
(void)W_prev;
(void)dt;
(void)gv;
(void)P_internal_data; // avoid warnings for unused params
}
// subtract all of the internal polarizations from the given f_minus_p
// field. Also given the fields array if it is needed for some reason.
// Only update for ft fields.
virtual void subtract_P(field_type ft, realnum *f_minus_p[NUM_FIELD_COMPONENTS][2],
void *P_internal_data) const {
(void)ft;
(void)f_minus_p;
(void)P_internal_data;
}
// whether, for the given field W, Meep needs to allocate P[c]
virtual bool needs_P(component c, int cmp, realnum *W[NUM_FIELD_COMPONENTS][2]) const;
// whether update_P will need the notowned part of W for this c
// (which means that Meep will need to communicate it between chunks)
virtual bool needs_W_notowned(component c, realnum *W[NUM_FIELD_COMPONENTS][2]) const;
// whether update_P needs the W_prev field (from the previous timestep)
virtual bool needs_W_prev() const { return false; }
/* A susceptibility may be associated with any amount of internal
data need to update the polarization field. This includes the
polarization field(s) itself. It may also, for example, store
the polarization field from previous timesteps, atomic-level
populations, or other data. These routines return the size of
this internal-data array and initialize it. */
virtual void *new_internal_data(realnum *W[NUM_FIELD_COMPONENTS][2],
const grid_volume &gv) const {
(void)W;
(void)gv;
return 0;
}
virtual void delete_internal_data(void *data) const;
virtual void init_internal_data(realnum *W[NUM_FIELD_COMPONENTS][2], realnum dt,
const grid_volume &gv, void *data) const {
(void)W;
(void)dt;
(void)gv;
(void)data;
}
virtual void *copy_internal_data(void *data) const {
(void)data;
return 0;
}
/* The following methods are used in boundaries.cpp to set up any
extra communications that may be necessary at chunk boundaries
for the internal data of a susceptibility's polarization
state. */
/* the number of notowned fields/data in the internal data that
are needed by update_P for the c Yee grid (note: we assume that we only
have internal data for c's where we have external polarizations) */
virtual int num_internal_notowned_needed(component c, void *P_internal_data) const {
(void)c;
(void)P_internal_data;
return 0;
}
/* the offset into the internal data of the n'th Yee-grid point in
the c Yee grid for the inotowned internal field, where
0 <= inotowned < size_internal_notowned_needed. */
virtual realnum *internal_notowned_ptr(int inotowned, component c, int n,
void *P_internal_data) const {
(void)inotowned;
(void)n;
(void)c;
(void)P_internal_data;
return 0;
}
/* same thing as above, except this gives (possibly complex)
internal fields that need to be multiplied by the same phase
factor as the fields at boundaries. Note: we assume internal fields
are complex if and only if !is_real (i.e. if EM fields are complex) */
virtual int num_cinternal_notowned_needed(component c, void *P_internal_data) const {
(void)c;
(void)P_internal_data;
return 0;
}
// real/imaginary parts offsets for cmp = 0/1
virtual realnum *cinternal_notowned_ptr(int inotowned, component c, int cmp, int n,
void *P_internal_data) const {
(void)inotowned;
(void)n;
(void)c;
(void)cmp;
(void)P_internal_data;
return 0;
}
virtual void dump_params(h5file *h5f, size_t *start) {
(void)h5f;
(void)start;
}
virtual int get_num_params() { return 0; }
// This should only be used when dumping and loading susceptibility data to hdf5
void set_id(int new_id) { id = new_id; };
susceptibility *next;
size_t ntot;
realnum *sigma[NUM_FIELD_COMPONENTS][5];
/* trivial_sigma[c][d] is true only if *none* of the processes has a
nontrivial sigma (c,d) component. This differs, from sigma,
which is non-NULL only if *this* process needs a nontrivial sigma
(c,d). Coordinated between processes at add_susceptibility, no
communication elsewhere. (We need this for boundary
communcations between chunks, where one chunk might have sigma ==
0 and the other != 0.) */
bool trivial_sigma[NUM_FIELD_COMPONENTS][5];
private:
static int cur_id; // unique id to assign to next susceptibility object
int id; // id for this object and its clones, for comparison purposes
};
/* a Lorentzian susceptibility
\chi(\omega) = sigma * omega_0^2 / (\omega_0^2 - \omega^2 - i\gamma \omega)
If no_omega_0_denominator is true, then we omit the omega_0^2 factor in the
denominator to obtain a Drude model. */
class lorentzian_susceptibility : public susceptibility {
public:
lorentzian_susceptibility(realnum omega_0, realnum gamma, bool no_omega_0_denominator = false)
: omega_0(omega_0), gamma(gamma), no_omega_0_denominator(no_omega_0_denominator) {}
virtual susceptibility *clone() const { return new lorentzian_susceptibility(*this); }
virtual ~lorentzian_susceptibility() {}
// Returns the 1st order nonlinear susceptibility
virtual std::complex<realnum> chi1(realnum freq, realnum sigma = 1);
virtual void update_P(realnum *W[NUM_FIELD_COMPONENTS][2],
realnum *W_prev[NUM_FIELD_COMPONENTS][2], realnum dt, const grid_volume &gv,
void *P_internal_data) const;
virtual void subtract_P(field_type ft, realnum *f_minus_p[NUM_FIELD_COMPONENTS][2],
void *P_internal_data) const;
virtual void *new_internal_data(realnum *W[NUM_FIELD_COMPONENTS][2], const grid_volume &gv) const;
virtual void init_internal_data(realnum *W[NUM_FIELD_COMPONENTS][2], realnum dt,
const grid_volume &gv, void *data) const;
virtual void *copy_internal_data(void *data) const;
virtual int num_cinternal_notowned_needed(component c, void *P_internal_data) const;
virtual realnum *cinternal_notowned_ptr(int inotowned, component c, int cmp, int n,
void *P_internal_data) const;
virtual void dump_params(h5file *h5f, size_t *start);
virtual int get_num_params() { return 4; }
protected:
realnum omega_0, gamma;
bool no_omega_0_denominator;
};
/* like a Lorentzian susceptibility, but the polarization equation
includes white noise with a specified amplitude */
class noisy_lorentzian_susceptibility : public lorentzian_susceptibility {
public:
noisy_lorentzian_susceptibility(realnum noise_amp, realnum omega_0, realnum gamma,
bool no_omega_0_denominator = false)
: lorentzian_susceptibility(omega_0, gamma, no_omega_0_denominator), noise_amp(noise_amp) {}
virtual susceptibility *clone() const { return new noisy_lorentzian_susceptibility(*this); }
virtual void update_P(realnum *W[NUM_FIELD_COMPONENTS][2],
realnum *W_prev[NUM_FIELD_COMPONENTS][2], realnum dt, const grid_volume &gv,
void *P_internal_data) const;
virtual void dump_params(h5file *h5f, size_t *start);
virtual int get_num_params() { return 5; }
protected:
realnum noise_amp;
};
typedef enum { GYROTROPIC_LORENTZIAN, GYROTROPIC_DRUDE, GYROTROPIC_SATURATED } gyrotropy_model;
/* gyrotropic susceptibility */
class gyrotropic_susceptibility : public susceptibility {
public:
gyrotropic_susceptibility(const vec &bias, realnum omega_0, realnum gamma, realnum alpha = 0.0,
gyrotropy_model model = GYROTROPIC_LORENTZIAN);
virtual susceptibility *clone() const { return new gyrotropic_susceptibility(*this); }
virtual void *new_internal_data(realnum *W[NUM_FIELD_COMPONENTS][2], const grid_volume &gv) const;
virtual void init_internal_data(realnum *W[NUM_FIELD_COMPONENTS][2], realnum dt,
const grid_volume &gv, void *data) const;
virtual void *copy_internal_data(void *data) const;
virtual bool needs_P(component c, int cmp, realnum *W[NUM_FIELD_COMPONENTS][2]) const;
virtual void update_P(realnum *W[NUM_FIELD_COMPONENTS][2],
realnum *W_prev[NUM_FIELD_COMPONENTS][2], realnum dt, const grid_volume &gv,
void *P_internal_data) const;
virtual void subtract_P(field_type ft, realnum *f_minus_p[NUM_FIELD_COMPONENTS][2],
void *P_internal_data) const;
virtual int num_cinternal_notowned_needed(component c, void *P_internal_data) const;
virtual realnum *cinternal_notowned_ptr(int inotowned, component c, int cmp, int n,
void *P_internal_data) const;
virtual void dump_params(h5file *h5f, size_t *start);
virtual int get_num_params() { return 8; }
virtual bool needs_W_notowned(component c, realnum *W[NUM_FIELD_COMPONENTS][2]) const {
(void)c;
(void)W;
return true;
}
protected:
realnum gyro_tensor[3][3];
realnum omega_0, gamma, alpha;
gyrotropy_model model;
};
class multilevel_susceptibility : public susceptibility {
public:
multilevel_susceptibility() : L(0), T(0), Gamma(0), N0(0), alpha(0), omega(0), gamma(0) {}
multilevel_susceptibility(int L, int T, const realnum *Gamma, const realnum *N0,
const realnum *alpha, const realnum *omega, const realnum *gamma,
const realnum *sigmat);
multilevel_susceptibility(const multilevel_susceptibility &from);
virtual susceptibility *clone() const { return new multilevel_susceptibility(*this); }
virtual ~multilevel_susceptibility();
virtual void update_P(realnum *W[NUM_FIELD_COMPONENTS][2],
realnum *W_prev[NUM_FIELD_COMPONENTS][2], realnum dt, const grid_volume &gv,
void *P_internal_data) const;
virtual void subtract_P(field_type ft, realnum *f_minus_p[NUM_FIELD_COMPONENTS][2],
void *P_internal_data) const;
virtual void *new_internal_data(realnum *W[NUM_FIELD_COMPONENTS][2], const grid_volume &gv) const;
virtual void init_internal_data(realnum *W[NUM_FIELD_COMPONENTS][2], realnum dt,
const grid_volume &gv, void *data) const;
virtual void *copy_internal_data(void *data) const;
virtual void delete_internal_data(void *data) const;
virtual int num_cinternal_notowned_needed(component c, void *P_internal_data) const;
virtual realnum *cinternal_notowned_ptr(int inotowned, component c, int cmp, int n,
void *P_internal_data) const;
// always need notowned W and W_prev for E dot dP/dt terms
virtual bool needs_W_notowned(component c, realnum *W[NUM_FIELD_COMPONENTS][2]) const {
(void)c;
(void)W;
return true;
}
virtual bool needs_W_prev() const { return true; }
protected:
int L; // number of atom levels
int T; // number of optical transitions
realnum *Gamma; // LxL matrix of relaxation rates Gamma[i*L+j] from i -> j
realnum *N0; // L initial populations
realnum *alpha; // LxT matrix of transition coefficients 1/omega
realnum *omega; // T transition frequencies
realnum *gamma; // T optical loss rates
realnum *sigmat; // 5*T transition-specific sigma-diagonal factors
};
// h5file.cpp: HDF5 file I/O. Most users, if they use this
// class at all, will only use the constructor to open the file, and
// will otherwise use the fields::output_hdf5 functions.
class h5file {
public:
typedef enum { READONLY, READWRITE, WRITE } access_mode;
h5file(const char *filename_, access_mode m = READWRITE, bool parallel_ = true);
~h5file(); // closes the files (and any open dataset)
bool ok();
void *read(const char *dataname, int *rank, size_t *dims, int maxrank,
bool single_precision = true);
void write(const char *dataname, int rank, const size_t *dims, void *data,
bool single_precision = true);
char *read(const char *dataname);
void write(const char *dataname, const char *data);
void create_data(const char *dataname, int rank, const size_t *dims, bool append_data = false,
bool single_precision = true);
void extend_data(const char *dataname, int rank, const size_t *dims);
void create_or_extend_data(const char *dataname, int rank, const size_t *dims, bool append_data,
bool single_precision = true);
void write_chunk(int rank, const size_t *chunk_start, const size_t *chunk_dims, float *data);
void write_chunk(int rank, const size_t *chunk_start, const size_t *chunk_dims, double *data);
void write_chunk(int rank, const size_t *chunk_start, const size_t *chunk_dims, size_t *data);
void done_writing_chunks();
void read_size(const char *dataname, int *rank, size_t *dims, int maxrank);
void read_chunk(int rank, const size_t *chunk_start, const size_t *chunk_dims, float *data);
void read_chunk(int rank, const size_t *chunk_start, const size_t *chunk_dims, double *data);
void read_chunk(int rank, const size_t *chunk_start, const size_t *chunk_dims, size_t *data);
void remove();
void remove_data(const char *dataname);
const char *file_name() const { return filename; }
void prevent_deadlock(); // hackery for exclusive mode
bool dataset_exists(const char *name);
private:
access_mode mode;
char *filename;
bool parallel;
bool is_cur(const char *dataname);
void unset_cur();
void set_cur(const char *dataname, void *data_id);
char *cur_dataname;
/* store hid_t values as hid_t* cast to void*, so that
files including meep.h don't need hdf5.h */
void *id; /* file */
void *cur_id; /* dataset, if any */
void *get_id(); // get current (file) id, opening/creating file if needed
void close_id();
public:
/* linked list to keep track of which datasets we are extending...
this is necessary so that create_or_extend_data can know whether
to create (overwrite) a dataset or extend it. */
struct extending_s {
int dindex;
char *dataname;
struct extending_s *next;
} * extending;
extending_s *get_extending(const char *dataname) const;
};
typedef double (*pml_profile_func)(double u, void *func_data);
#define DEFAULT_SUBPIXEL_TOL 1e-4
#define DEFAULT_SUBPIXEL_MAXEVAL 100000
/* This class is used to compute position-dependent material properties
like the dielectric function, permeability (mu), polarizability sigma,
nonlinearities, et cetera. Simple cases of stateless functions are
handled by canned subclasses below, but more complicated cases
can be handled by creating a user-defined subclass of material_function.
It is useful to group different properties into one class because
it is likely that complicated implementations will share state between
properties. */
class material_function {
material_function(const material_function &ef) { (void)ef; } // prevent copying
public:
material_function() {}
virtual ~material_function() {}
/* Specify a restricted grid_volume: all subsequent eps/sigma/etc
calls will be for points inside v, until the next set_volume. */
virtual void set_volume(const volume &v) { (void)v; }
virtual void unset_volume(void) {} // unrestrict the grid_volume
virtual double chi1p1(field_type ft, const vec &r) {
(void)ft;
(void)r;
return 1.0;
}
/* scalar dielectric function */
virtual double eps(const vec &r) { return chi1p1(E_stuff, r); }
/* scalar permeability function */
virtual bool has_mu() { return false; } /* true if mu != 1 */
virtual double mu(const vec &r) { return chi1p1(H_stuff, r); }
/* scalar conductivity function */
virtual bool has_conductivity(component c) {
(void)c;
return false;
}
virtual double conductivity(component c, const vec &r) {
(void)c;
(void)r;
return 0.0;
}
// fallback routine based on spherical quadrature
vec normal_vector(field_type ft, const volume &v);
/* Return c'th row of effective 1/(1+chi1) tensor in the given grid_volume v
... virtual so that e.g. libctl can override with more-efficient
libctlgeom-based routines. maxeval == 0 if no averaging desired. */
virtual void eff_chi1inv_row(component c, double chi1inv_row[3], const volume &v,
double tol = DEFAULT_SUBPIXEL_TOL,
int maxeval = DEFAULT_SUBPIXEL_MAXEVAL);
/* polarizability sigma function: return c'th row of tensor */
virtual void sigma_row(component c, double sigrow[3], const vec &r) {
(void)c;
(void)r;
sigrow[0] = sigrow[1] = sigrow[2] = 0.0;
}
// Nonlinear susceptibilities
virtual bool has_chi3(component c) {
(void)c;
return false;
}
virtual double chi3(component c, const vec &r) {
(void)c;
(void)r;
return 0.0;
}
virtual bool has_chi2(component c) {
(void)c;
return false;
}
virtual double chi2(component c, const vec &r) {
(void)c;
(void)r;
return 0.0;
}
};
class simple_material_function : public material_function {
double (*f)(const vec &);
public:
simple_material_function(double (*func)(const vec &)) { f = func; }
virtual ~simple_material_function() {}
virtual double chi1p1(field_type ft, const vec &r) {
(void)ft;
return f(r);
}
virtual double eps(const vec &r) { return f(r); }
virtual double mu(const vec &r) { return f(r); }
virtual double conductivity(component c, const vec &r) {
(void)c;
return f(r);
}
virtual void sigma_row(component c, double sigrow[3], const vec &r) {
sigrow[0] = sigrow[1] = sigrow[2] = 0.0;
sigrow[component_index(c)] = f(r);
}
virtual double chi3(component c, const vec &r) {
(void)c;
return f(r);
}
virtual double chi2(component c, const vec &r) {
(void)c;
return f(r);
}
};
class structure;
class structure_chunk {
public:
double a, Courant, dt; // resolution a, Courant number, and timestep dt=Courant/a
realnum *chi3[NUM_FIELD_COMPONENTS], *chi2[NUM_FIELD_COMPONENTS];
realnum *chi1inv[NUM_FIELD_COMPONENTS][5];
bool trivial_chi1inv[NUM_FIELD_COMPONENTS][5];
realnum *conductivity[NUM_FIELD_COMPONENTS][5];
realnum *condinv[NUM_FIELD_COMPONENTS][5]; // cache of 1/(1+conduct*dt/2)
bool condinv_stale; // true if condinv needs to be recomputed
realnum *sig[6], *kap[6], *siginv[6]; // conductivity array for uPML
int sigsize[6]; // conductivity array size
grid_volume gv; // integer grid_volume that could be bigger than non-overlapping v below
volume v;
susceptibility *chiP[NUM_FIELD_TYPES]; // only E_stuff and H_stuff are used
double
cost; // The cost of this chunk's grid_volume as computed by split_by_cost and fragment_stats
int refcount; // reference count of objects using this structure_chunk
~structure_chunk();
structure_chunk(const grid_volume &gv, const volume &vol_limit, double Courant, int proc_num);
structure_chunk(const structure_chunk *);
void set_chi1inv(component c, material_function &eps, bool use_anisotropic_averaging, double tol,
int maxeval);
bool has_chi(component c, direction d) const;
bool has_chisigma(component c, direction d) const;
bool has_chi1inv(component c, direction d) const;
void set_conductivity(component c, material_function &eps);
void update_condinv();
void set_chi3(component c, material_function &eps);
void set_chi2(component c, material_function &eps);
void use_pml(direction, double dx, double boundary_loc, double Rasymptotic, double mean_stretch,
pml_profile_func pml_profile, void *pml_profile_data, double pml_profile_integral,
double pml_profile_integral_u);
void add_susceptibility(material_function &sigma, field_type ft, const susceptibility &sus);
void mix_with(const structure_chunk *, double);
int n_proc() const { return the_proc; } // Says which proc owns me!
int is_mine() const { return the_is_mine; }
void remove_susceptibilities();
// monitor.cpp
std::complex<double> get_chi1inv_at_pt(component, direction, int idx, double frequency = 0) const;
std::complex<double> get_chi1inv(component, direction, const ivec &iloc,
double frequency = 0) const;
std::complex<double> get_inveps(component c, direction d, const ivec &iloc,
double frequency = 0) const {
return get_chi1inv(c, d, iloc, frequency);
}
double max_eps() const;
private:
double pml_fmin;
int the_proc;
int the_is_mine;
};
// linked list of descriptors for boundary regions (currently just for PML)
class boundary_region {
public:
typedef enum { NOTHING_SPECIAL, PML } boundary_region_kind;
boundary_region()
: kind(NOTHING_SPECIAL), thickness(0.0), Rasymptotic(1e-16), mean_stretch(1.0),
pml_profile(NULL), pml_profile_data(NULL), pml_profile_integral(1.0),
pml_profile_integral_u(1.0), d(NO_DIRECTION), side(Low), next(0) {}
boundary_region(boundary_region_kind kind, double thickness, double Rasymptotic,
double mean_stretch, pml_profile_func pml_profile, void *pml_profile_data,
double pml_profile_integral, double pml_profile_integral_u, direction d,
boundary_side side, boundary_region *next = 0)
: kind(kind), thickness(thickness), Rasymptotic(Rasymptotic), mean_stretch(mean_stretch),
pml_profile(pml_profile), pml_profile_data(pml_profile_data),
pml_profile_integral(pml_profile_integral), pml_profile_integral_u(pml_profile_integral_u),
d(d), side(side), next(next) {}
boundary_region(const boundary_region &r)
: kind(r.kind), thickness(r.thickness), Rasymptotic(r.Rasymptotic),
mean_stretch(r.mean_stretch), pml_profile(r.pml_profile),
pml_profile_data(r.pml_profile_data), pml_profile_integral(r.pml_profile_integral),
pml_profile_integral_u(r.pml_profile_integral_u), d(r.d), side(r.side) {
next = r.next ? new boundary_region(*r.next) : 0;
}
~boundary_region() {
if (next) delete next;
}
void operator=(const boundary_region &r) {
kind = r.kind;
thickness = r.thickness;
Rasymptotic = r.Rasymptotic;
mean_stretch = r.mean_stretch;
pml_profile = r.pml_profile;
pml_profile_data = r.pml_profile_data;
pml_profile_integral = r.pml_profile_integral;
pml_profile_integral_u = r.pml_profile_integral_u;
d = r.d;
side = r.side;
if (next) delete next;
next = r.next ? new boundary_region(*r.next) : 0;
}
boundary_region operator+(const boundary_region &r0) const {
boundary_region r(*this), *cur = &r;
while (cur->next)
cur = cur->next;
cur->next = new boundary_region(r0);
return r;
}
boundary_region operator*(double strength_mult) const {
boundary_region r(*this), *cur = &r;
while (cur) {
cur->Rasymptotic = pow(cur->Rasymptotic, strength_mult);
cur = cur->next;
}
return r;
}
void apply(structure *s) const;
void apply(const structure *s, structure_chunk *sc) const;
bool check_ok(const grid_volume &gv) const;
private:
boundary_region_kind kind;
double thickness, Rasymptotic, mean_stretch;
pml_profile_func pml_profile;
void *pml_profile_data;
double pml_profile_integral, pml_profile_integral_u;
direction d;
boundary_side side;
boundary_region *next;
};
boundary_region pml(double thickness, direction d, boundary_side side, double Rasymptotic = 1e-15,
double mean_stretch = 1.0);
boundary_region pml(double thickness, direction d, double Rasymptotic = 1e-15,
double mean_stretch = 1.0);
boundary_region pml(double thickness, double Rasymptotic = 1e-15, double mean_stretch = 1.0);
#define no_pml() boundary_region()
class binary_partition;
enum in_or_out { Incoming = 0, Outgoing };
const std::initializer_list<in_or_out> all_in_or_out{Incoming, Outgoing};
enum connect_phase { CONNECT_PHASE = 0, CONNECT_NEGATE = 1, CONNECT_COPY = 2 };
const std::initializer_list<connect_phase> all_connect_phases{
CONNECT_PHASE, CONNECT_NEGATE, CONNECT_COPY};
constexpr int NUM_CONNECT_PHASE_TYPES = 3;
// Communication pair(i, j) implies data being sent from i to j.
using chunk_pair = std::pair<int, int>;
// Key for the map that stored communication sizes.
struct comms_key {
field_type ft;
connect_phase phase;
chunk_pair pair;
};
// Defined in fields.cpp
bool operator==(const comms_key &lhs, const comms_key &rhs);
class comms_key_hash_fn {
public:
inline std::size_t operator()(const comms_key &key) const {
// Unroll hash combination to promote the generatiion of efficient code.
std::size_t ret = int_hasher(int(key.ft));
ret ^= int_hasher(int(key.phase)) + kHashAddConst + (ret << 6) + (ret >> 2);
ret ^= int_hasher(key.pair.first) + kHashAddConst + (ret << 6) + (ret >> 2);
ret ^= int_hasher(key.pair.second) + kHashAddConst + (ret << 6) + (ret >> 2);
return ret;
}
private:
static constexpr size_t kHashAddConst = 0x9e3779b9;
std::hash<int> int_hasher;
};
// Represents a communication operation between chunks.
struct comms_operation {
ptrdiff_t my_chunk_idx;
ptrdiff_t other_chunk_idx;
int other_proc_id;
int pair_idx; // The numeric pair index used in many communications related arrays.
size_t transfer_size;
in_or_out comm_direction;
int tag;
};
struct comms_sequence {
std::vector<comms_operation> receive_ops;
std::vector<comms_operation> send_ops;
void clear() {
receive_ops.clear();
send_ops.clear();
}
};
// RAII based comms_manager that allows asynchronous send and receive functions to be initiated.
// Upon destruction, the comms_manager waits for completion of all enqueued operations.
class comms_manager {
public:
virtual ~comms_manager() {}
virtual void send_real_async(const void *buf, size_t count, int dest, int tag) = 0;
virtual void receive_real_async(void *buf, size_t count, int source, int tag) = 0;
virtual size_t max_transfer_size() const { return std::numeric_limits<size_t>::max(); };
};
// Factory function for `comms_manager`.
std::unique_ptr<comms_manager> create_comms_manager();
class structure {
public:
structure_chunk **chunks;
int num_chunks;
bool shared_chunks; // whether modifications to chunks will be visible to fields objects
grid_volume gv, user_volume;
double a, Courant, dt; // resolution a, Courant number, and timestep dt=Courant/a
volume v;
symmetry S;
const char *outdir;
grid_volume *effort_volumes;
double *effort;
int num_effort_volumes;
~structure();
structure(const grid_volume &gv, material_function &eps,
const boundary_region &br = boundary_region(), const symmetry &s = meep::identity(),
int num_chunks = 0, double Courant = 0.5, bool use_anisotropic_averaging = false,
double tol = DEFAULT_SUBPIXEL_TOL, int maxeval = DEFAULT_SUBPIXEL_MAXEVAL,
const binary_partition *_bp = NULL);
structure(const grid_volume &gv, double eps(const vec &),
const boundary_region &br = boundary_region(), const symmetry &s = meep::identity(),
int num_chunks = 0, double Courant = 0.5, bool use_anisotropic_averaging = false,
double tol = DEFAULT_SUBPIXEL_TOL, int maxeval = DEFAULT_SUBPIXEL_MAXEVAL,
const binary_partition *_bp = NULL);
structure(const structure &);
void set_materials(material_function &mat, bool use_anisotropic_averaging = true,
double tol = DEFAULT_SUBPIXEL_TOL, int maxeval = DEFAULT_SUBPIXEL_MAXEVAL);
void set_chi1inv(component c, material_function &eps, bool use_anisotropic_averaging = true,
double tol = DEFAULT_SUBPIXEL_TOL, int maxeval = DEFAULT_SUBPIXEL_MAXEVAL);
bool has_chi(component c, direction d) const;
void set_epsilon(material_function &eps, bool use_anisotropic_averaging = true,
double tol = DEFAULT_SUBPIXEL_TOL, int maxeval = DEFAULT_SUBPIXEL_MAXEVAL);
void set_epsilon(double eps(const vec &), bool use_anisotropic_averaging = true,
double tol = DEFAULT_SUBPIXEL_TOL, int maxeval = DEFAULT_SUBPIXEL_MAXEVAL);
void set_mu(material_function &eps, bool use_anisotropic_averaging = true,
double tol = DEFAULT_SUBPIXEL_TOL, int maxeval = DEFAULT_SUBPIXEL_MAXEVAL);
void set_mu(double mu(const vec &), bool use_anisotropic_averaging = true,
double tol = DEFAULT_SUBPIXEL_TOL, int maxeval = DEFAULT_SUBPIXEL_MAXEVAL);
void set_conductivity(component c, material_function &conductivity);
void set_conductivity(component C, double conductivity(const vec &));
void set_chi3(component c, material_function &eps);
void set_chi3(material_function &eps);
void set_chi3(double eps(const vec &));
void set_chi2(component c, material_function &eps);
void set_chi2(material_function &eps);
void set_chi2(double eps(const vec &));
void add_susceptibility(double sigma(const vec &), field_type c, const susceptibility &sus);
void add_susceptibility(material_function &sigma, field_type c, const susceptibility &sus);
void remove_susceptibilities();
void set_output_directory(const char *name);
void mix_with(const structure *, double);
bool equal_layout(const structure &) const;
void print_layout(void) const;
std::vector<grid_volume> get_chunk_volumes() const;
std::vector<int> get_chunk_owners() const;
// structure_dump.cpp
void dump(const char *filename);
void dump_chunk_layout(const char *filename);
void load(const char *filename);
void load_chunk_layout(const char *filename, boundary_region &br);
void load_chunk_layout(const std::vector<grid_volume> &gvs,
const std::vector<int> &ids,
boundary_region &br);
// monitor.cpp
std::complex<double> get_chi1inv(component, direction, const ivec &origloc, double frequency = 0,
bool parallel = true) const;
std::complex<double> get_chi1inv(component, direction, const vec &loc, double frequency = 0,
bool parallel = true) const;
std::complex<double> get_inveps(component c, direction d, const ivec &origloc,
double frequency = 0) const {
return get_chi1inv(c, d, origloc, frequency);
}
std::complex<double> get_inveps(component c, direction d, const vec &loc,
double frequency = 0) const {
return get_chi1inv(c, d, loc, frequency);
}
std::complex<double> get_eps(const vec &loc, double frequency = 0) const;
std::complex<double> get_mu(const vec &loc, double frequency = 0) const;
double max_eps() const;
double estimated_cost(int process = my_rank());
// Returns the binary partition that was used to partition the volume into chunks. The returned
// pointer is only valid for the lifetime of this `structure` instance.
const binary_partition *get_binary_partition() const;
friend class boundary_region;
private:
void use_pml(direction d, boundary_side b, double dx);
void add_to_effort_volumes(const grid_volume &new_effort_volume, double extra_effort);
void choose_chunkdivision(const grid_volume &gv, int num_chunks, const boundary_region &br,
const symmetry &s, const binary_partition *_bp);
void check_chunks();
void changing_chunks();
// Helper methods for dumping and loading susceptibilities
void set_chiP_from_file(h5file *file, const char *dataset, field_type ft);
void write_susceptibility_params(h5file *file, const char *dname, int EorH);
std::unique_ptr<binary_partition> bp;
};
// defined in structure.cpp
std::unique_ptr<binary_partition> choose_chunkdivision(grid_volume &gv, volume &v, int num_chunks,
const symmetry &s);
// defined in structure_dump.cpp
void split_by_binarytree(grid_volume gvol,
std::vector<grid_volume> &result_gvs,
std::vector<int> &result_ids,
const binary_partition *bp);
class src_vol;
class fields;
class fields_chunk;
class flux_vol;
// Time-dependence of a current source, intended to be overridden by
// subclasses. current() and dipole() are be related by
// current = d(dipole)/dt (or rather, the finite-difference equivalent).
class src_time {
public:
// the following variable specifies whether the current
// source is specified as a current or as an integrated
// current (a dipole moment), if possible. In the original Meep,
// by default electric sources are integrated and magnetic
// sources are not, but this may change.
bool is_integrated;
src_time() {
is_integrated = true;
current_time = nan;
current_current = 0.0;
next = NULL;
}
virtual ~src_time() { delete next; }
src_time(const src_time &t) {
is_integrated = t.is_integrated;
current_time = t.current_time;
current_current = t.current_current;
current_dipole = t.current_dipole;
if (t.next)
next = t.next->clone();
else
next = NULL;
}
std::complex<double> dipole() const { return current_dipole; }
std::complex<double> current() const { return current_current; }
void update(double time, double dt) {
if (time != current_time) {
current_dipole = dipole(time);
current_current = current(time, dt);
current_time = time;
}
}
// subclasses *can* override this method in order to specify the
// current directly rather than as the derivative of dipole.
// in that case you would probably ignore the dt argument.
virtual std::complex<double> current(double time, double dt) const {
return ((dipole(time + dt) - dipole(time)) / dt);
}
double last_time_max() { return last_time_max(0.0); }
double last_time_max(double after);
src_time *add_to(src_time *others, src_time **added) const;
src_time *next;
// subclasses should override these methods:
virtual std::complex<double> dipole(double time) const {
(void)time;
return 0;
}
virtual double last_time() const { return 0.0; }
virtual src_time *clone() const { return new src_time(*this); }
virtual bool is_equal(const src_time &t) const {
(void)t;
return 1;
}
virtual std::complex<double> frequency() const { return 0.0; }
virtual void set_frequency(std::complex<double> f) { (void)f; }
private:
double current_time;
std::complex<double> current_dipole, current_current;
};
bool src_times_equal(const src_time &t1, const src_time &t2);
// Gaussian-envelope source with given frequency, width, peak-time, cutoff
class gaussian_src_time : public src_time {
public:
gaussian_src_time(double f, double fwidth, double s = 5.0);
gaussian_src_time(double f, double w, double start_time, double end_time);
virtual ~gaussian_src_time() {}
virtual std::complex<double> dipole(double time) const;
virtual double last_time() const { return float(peak_time + cutoff); };
virtual src_time *clone() const { return new gaussian_src_time(*this); }
virtual bool is_equal(const src_time &t) const;
virtual std::complex<double> frequency() const { return freq; }
virtual void set_frequency(std::complex<double> f) { freq = real(f); }
std::complex<double> fourier_transform(const double f);
private: