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
|
#include <meep.hpp>
using namespace meep;
#include "config.h"
using namespace std;
const double diameter = 0.8;
const double r = diameter*0.5;
double holey_2d(const vec &xx) {
const grid_volume gv = vol2d(2.0,1.0,100.0);
vec p = xx - gv.center();
while (p.x() <-0.5) p += vec(1.0,0);
while (p.x() > 0.5) p -= vec(1.0,0);
while (p.y() <-0.5) p += vec(0,1.0);
while (p.y() > 0.5) p -= vec(0,1.0);
if (fabs(p & p) < r*r - 1e-12) return 1.0;
return 12.0;
}
double holey_shifted_2d(const vec &xx) {
return holey_2d(xx + vec(pi*0.01, 3 - pi)*0.5);
}
double get_the_freq(monitor_point *p, component c) {
complex<double> *amp, *freqs;
int num;
p->harminv(c, &, &freqs, &num, 0.15, 0.20, 8);
if (!num) return 0.0;
double best_amp = abs(amp[0]), best_freq = fabs(real(freqs[0]));
for (int i=1;i<num;i++)
if (abs(amp[i]) > best_amp &&
fabs(imag(freqs[i])/real(freqs[i])) < 0.002) {
best_amp = abs(amp[i]);
best_freq = fabs(real(freqs[i]));
}
delete[] freqs;
delete[] amp;
return best_freq;
}
double freq_at_resolution(double e(const vec &), double a, component c,
double beta) {
const grid_volume gv = vol2d(2.0,1.0,a);
structure s(gv, e);
s.set_epsilon(e);
fields f(&s, 0, beta);
f.use_real_fields();
f.use_bloch(vec(0,0));
f.add_point_source(c, 0.18, 2.5, 0.0, 6.0, vec(0.5,0.5), 1.0);
f.add_point_source(c, 0.18, 2.5, 0.0, 6.0, vec(1.5,0.5),-1.0);
while (f.time() <= f.last_source_time() + 10.0 && !interrupt) f.step();
const double fourier_timesteps = 3000.0;
const double ttot = fourier_timesteps/a + f.time();
monitor_point *p = NULL;
while (f.time() <= ttot) {
f.step();
p = f.get_new_point(vec(0.52,0.97), p);
}
const double freq = get_the_freq(p, c);
delete p;
return freq;
}
void check_convergence(component c, double best_guess, double beta)
{
const double amin = 5.0, amax = 30.0, adelta = 5.0;
master_printf("Checking convergence for %s field...\n",
component_name(c));
if (beta != 0)
master_printf("... using exp(i beta z) z-dependence with beta=%g\n", beta);
if (best_guess)
master_printf("(The correct frequency should be %g.)\n", best_guess);
for (double a=amax; a >= amin; a-=adelta) {
const double freq = freq_at_resolution(holey_2d, a, c, beta);
const double freq_shifted = freq_at_resolution(holey_shifted_2d, a, c,
beta);
// Initialize best guess at the correct freq.
if (!best_guess) {
best_guess = freq + 0.5*(freq_shifted - freq);
master_printf("The frequency is approximately %g\n", best_guess);
} else {
master_printf("frequency for a=%g is %g, %g (shifted), %g (mean)\n",
a, freq, freq_shifted, 0.5 * (freq + freq_shifted));
master_printf("Unshifted freq error is %g/%g/%g\n",
(freq - best_guess)*a*a, a, a);
if (fabs(freq - best_guess)*a*a > 0.4)
abort("Frequency doesn't converge properly with a.\n");
master_printf("Shifted freq error is %g/%g/%g\n",
(freq_shifted - best_guess)*a*a, a, a);
if (fabs(freq_shifted - best_guess)*a*a > 0.4)
abort("Frequency doesn't converge properly with a.\n");
}
// Check frequency difference...
master_printf("Frequency difference with a of %g is %g/%g/%g\n",
a, (freq - freq_shifted)*a*a, a, a);
if (fabs(freq - freq_shifted)*a*a > 0.4)
abort("Frequency difference = doesn't converge properly with a.\n");
}
master_printf("Passed 2D resolution convergence test for %s!\n",
component_name(c));
}
int main(int argc, char **argv) {
initialize mpi(argc, argv);
quiet = true;
#ifdef HAVE_HARMINV
master_printf("Running holes square-lattice resolution convergence test.\n");
check_convergence(Ey, 0.179944, 0); // from MPB; correct to >= 4 dec. places
check_convergence(Ez, 0.166998, 0); // from MPB; correct to >= 4 dec. places
check_convergence(Ez, 0.173605, .1); // from MPB; correct to >= 4 dec. places
#endif
return 0;
}
|
