/
bend-flux-ll.cpp
300 lines (261 loc) · 11 KB
/
bend-flux-ll.cpp
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
/***************************************************************/
/* C++ port of meep/examples/bend-flux.ctl, using the */
/* "low-level" meep C++ interface stack, which consists of */
/* libmeep_geom + libctlgeom + libmeep */
/***************************************************************/
/*
; From the Meep tutorial: transmission around a 90-degree waveguide
; bend in 2d.
*/
#include <stdio.h>
#include <stdlib.h>
#include <complex>
#include "meep.hpp"
#include "ctl-math.h"
#include "ctlgeom.h"
#include "meepgeom.hpp"
using namespace meep;
typedef std::complex<double> cdouble;
vector3 v3(double x, double y = 0.0, double z = 0.0) {
vector3 v;
v.x = x;
v.y = y;
v.z = z;
return v;
}
/***************************************************************/
/* dummy material function needed to pass to structure( ) */
/* constructor as a placeholder before we can call */
/* set_materials_from_geometry */
/***************************************************************/
double dummy_eps(const vec &) { return 1.0; }
/***************************************************************/
/* define geometry from GDSII file *****************************/
/***************************************************************/
#define GEOM_LAYER 0 // hard-coded indices of GDSII layers
#define STRAIGHT_WVG_LAYER 1 // on which various geometric entities are defined
#define BENT_WVG_LAYER 2
#define SOURCE_LAYER 3
#define RFLUX_LAYER 4
#define STRAIGHT_TFLUX_LAYER 5
#define BENT_TFLUX_LAYER 6
structure create_structure_from_GDSII(char *GDSIIFile, bool no_bend, volume &vsrc, volume &vrefl,
volume &vtrans) {
// set computational cell
double dpml = 1.0;
double resolution = 10.0;
grid_volume gv = meep_geom::set_geometry_from_GDSII(resolution, GDSIIFile, GEOM_LAYER);
structure the_structure(gv, dummy_eps, pml(dpml));
// define waveguide
geometric_object objects[1];
meep_geom::material_type dielectric = meep_geom::make_dielectric(12.0);
int GDSIILayer = (no_bend ? STRAIGHT_WVG_LAYER : BENT_WVG_LAYER);
objects[0] = meep_geom::get_GDSII_prism(dielectric, GDSIIFile, GDSIILayer);
geometric_object_list g = {1, objects};
meep_geom::set_materials_from_geometry(&the_structure, g);
// define volumes for source and flux-monitor regions
vsrc = meep_geom::get_GDSII_volume(GDSIIFile, SOURCE_LAYER);
vrefl = meep_geom::get_GDSII_volume(GDSIIFile, RFLUX_LAYER);
vtrans =
meep_geom::get_GDSII_volume(GDSIIFile, (no_bend ? STRAIGHT_TFLUX_LAYER : BENT_TFLUX_LAYER));
return the_structure;
}
/***************************************************************/
/***************************************************************/
/***************************************************************/
structure create_structure_by_hand(bool no_bend, bool use_prisms, volume &vsrc, volume &vrefl,
volume &vtrans) {
double sx = 16.0; // size of cell in X direction
double sy = 32.0; // size of cell in Y direction
double pad = 4.0; // padding distance between waveguide and cell edge
double w = 1.0; // width of waveguide
double dpml = 1.0; // PML thickness
double resolution = 10.0;
// (set! geometry-lattice (make lattice (size sx sy no-size)))
// (set! pml-layers (list (make pml (thickness 1.0))))
geometry_lattice.size.x = sx;
geometry_lattice.size.y = sy;
geometry_lattice.size.z = 0.0;
grid_volume gv = voltwo(sx, sy, resolution);
gv.center_origin();
structure the_structure(gv, dummy_eps, pml(dpml));
double wvg_ycen = -0.5 * (sy - w - 2.0 * pad); // y center of horiz. wvg
double wvg_xcen = 0.5 * (sx - w - 2.0 * pad); // x center of vert. wvg
// (set! geometry
// (if no-bend?
// (list
// (list
// (make block
// (center wvg-xcen (* 0.5 pad))
// (size w (- sy pad) infinity)
// (material (make dielectric (epsilon 12)))))))
// (make block
// (center 0 wvg-ycen)
// (size infinity w infinity)
// (material (make dielectric (epsilon 12)))))
// (make block
// (center (* -0.5 pad) wvg-ycen)
// (size (- sx pad) w infinity)
// (material (make dielectric (epsilon 12))))
vector3 e1 = v3(1.0, 0.0, 0.0);
vector3 e2 = v3(0.0, 1.0, 0.0);
vector3 e3 = v3(0.0, 0.0, 1.0);
meep_geom::material_type dielectric = meep_geom::make_dielectric(12.0);
if (no_bend) {
if (use_prisms) {
vector3 vertices[4];
vertices[0] = v3(-0.5 * sx, wvg_ycen - 0.5 * w, 0.0);
vertices[1] = v3(+0.5 * sx, wvg_ycen - 0.5 * w, 0.0);
vertices[2] = v3(+0.5 * sx, wvg_ycen + 0.5 * w, 0.0);
vertices[3] = v3(-0.5 * sx, wvg_ycen + 0.5 * w, 0.0);
double height = 0.0;
vector3 axis = v3(0.0, 0.0, 1.0);
geometric_object objects[1];
objects[0] = make_prism(dielectric, vertices, 4, height, axis);
geometric_object_list g = {1, objects};
meep_geom::set_materials_from_geometry(&the_structure, g);
}
else {
vector3 center = v3(0.0, wvg_ycen, 0.0);
vector3 size = v3(sx, w, 0.0);
geometric_object objects[1];
objects[0] = make_block(dielectric, center, e1, e2, e3, size);
geometric_object_list g = {1, objects};
meep_geom::set_materials_from_geometry(&the_structure, g);
}
}
else {
if (use_prisms) {
vector3 vertices[6];
vertices[0] = v3(-0.5 * sx, wvg_ycen - 0.5 * w);
vertices[1] = v3(wvg_xcen + 0.5 * w, wvg_ycen - 0.5 * w);
vertices[2] = v3(wvg_xcen + 0.5 * w, +0.5 * sy);
vertices[3] = v3(wvg_xcen - 0.5 * w, +0.5 * sy);
vertices[4] = v3(wvg_xcen - 0.5 * w, wvg_ycen + 0.5 * w);
vertices[5] = v3(-0.5 * sx, wvg_ycen + 0.5 * w);
double height = 0.0;
vector3 axis = v3(0.0, 0.0, 1.0);
geometric_object objects[1];
objects[0] = make_prism(dielectric, vertices, 6, height, axis);
geometric_object_list g = {1, objects};
meep_geom::set_materials_from_geometry(&the_structure, g);
}
else {
geometric_object objects[2];
vector3 hcenter = v3(-0.5 * pad, wvg_ycen), hsize = v3(sx - pad, w);
vector3 vcenter = v3(wvg_xcen, 0.5 * pad), vsize = v3(w, sy - pad);
objects[0] = make_block(dielectric, hcenter, e1, e2, e3, hsize);
objects[1] = make_block(dielectric, vcenter, e1, e2, e3, vsize);
geometric_object_list g = {2, objects};
meep_geom::set_materials_from_geometry(&the_structure, g);
}
}
vsrc =
volume(vec(-0.5 * sx + dpml, wvg_ycen - 0.5 * w), vec(-0.5 * sx + dpml, wvg_ycen + 0.5 * w));
vrefl = volume(vec(-0.5 * sx + 1.5, wvg_ycen - w), vec(-0.5 * sx + 1.5, wvg_ycen + w));
vtrans = (no_bend ? volume(vec(0.5 * sx - 1.5, wvg_ycen - w), vec(0.5 * sx - 1.5, wvg_ycen + w))
: volume(vec(wvg_xcen - w, 0.5 * sy - 1.5), vec(wvg_xcen + w, 0.5 * sy - 1.5)));
return the_structure;
}
/***************************************************************/
/***************************************************************/
/***************************************************************/
void bend_flux(bool no_bend, char *GDSIIFile, bool use_prisms) {
vec v0;
volume vsrc(v0), vrefl(v0), vtrans(v0);
structure the_structure =
GDSIIFile ? create_structure_from_GDSII(GDSIIFile, no_bend, vsrc, vrefl, vtrans)
: create_structure_by_hand(no_bend, use_prisms, vsrc, vrefl, vtrans);
fields f(&the_structure);
// (set! sources (list
// (make source
// (src (make gaussian-src (frequency fcen) (fwidth df)))
// (component Ez)
// (center (+ 1 (* -0.5 sx)) wvg-ycen)
// (size 0 w))))
//
double fcen = 0.15; // ; pulse center frequency
double df = 0.1; // ; df
gaussian_src_time src(fcen, df);
f.add_volume_source(Ez, src, vsrc);
//(define trans ; transmitted flux
// (add-flux fcen df nfreq
// (if no-bend?
// (make flux-region
// (center (- (/ sx 2) 1.5) wvg-ycen) (size 0 (* w 2)))
// (make flux-region
// (center wvg-xcen (- (/ sy 2) 1.5)) (size (* w 2) 0)))))
double f_start = fcen - 0.5 * df;
double f_end = fcen + 0.5 * df;
int nfreq = 100; // number of frequencies at which to compute flux
direction trans_dir = no_bend ? X : Y;
dft_flux trans = f.add_dft_flux(trans_dir, vtrans, f_start, f_end, nfreq);
//(define refl ; reflected flux
// (add-flux fcen df nfreq
// (make flux-region
// (center (+ (* -0.5 sx) 1.5) wvg-ycen) (size 0 (* w 2)))))
//
dft_flux refl = f.add_dft_flux(NO_DIRECTION, vrefl, f_start, f_end, nfreq);
char *prefix = const_cast<char *>(no_bend ? "straight" : "bent");
f.output_hdf5(Dielectric, f.total_volume(), 0, false, true, prefix);
//; for normal run, load negated fields to subtract incident from refl. fields
//(if (not no-bend?) (load-minus-flux "refl-flux" refl))
const char *dataname = "refl-flux";
if (!no_bend) {
refl.load_hdf5(f, dataname);
refl.scale_dfts(-1.0);
}
//(run-sources+
// (stop-when-fields-decayed 50 Ez
// (if no-bend?
// (vector3 (- (/ sx 2) 1.5) wvg-ycen)
// (vector3 wvg-xcen (- (/ sy 2) 1.5)))
// 1e-3)
// (at-beginning output-epsilon))
vec eval_point = vtrans.center();
double DeltaT = 10.0, NextCheckTime = f.round_time() + DeltaT;
double Tol = 1.0e-3;
double max_abs = 0.0, cur_max = 0.0;
bool Done = false;
do {
f.step();
// manually check fields-decayed condition
double absEz = abs(f.get_field(Ez, eval_point));
cur_max = fmax(cur_max, absEz);
if (f.round_time() >= NextCheckTime) {
NextCheckTime += DeltaT;
max_abs = fmax(max_abs, cur_max);
if ((max_abs > 0.0) && cur_max < Tol * max_abs) Done = true;
cur_max = 0.0;
}
// printf("%.2e %.2e %.2e %.2e\n",f.round_time(),absEz,max_abs,cur_max);
} while (!Done);
//; for normalization run, save flux fields for refl. plane
//(if no-bend? (save-flux "refl-flux" refl))
if (no_bend) refl.save_hdf5(f, dataname);
//(display-fluxes trans refl)
printf("%11s | %12s | %12s\n", " Freq ", " trans flux", " refl flux");
double f0 = fcen - 0.5 * df, fstep = df / (nfreq - 1);
double *trans_flux = trans.flux();
double *refl_flux = refl.flux();
for (int nf = 0; nf < nfreq; nf++)
printf("%.4e | %+.4e | %+.4e\n", f0 + nf * fstep, trans_flux[nf], refl_flux[nf]);
}
/***************************************************************/
/***************************************************************/
/***************************************************************/
int main(int argc, char *argv[]) {
initialize mpi(argc, argv);
bool use_prisms = false;
char *GDSIIFile = 0;
for (int narg = 1; narg < argc; narg++)
if (!strcasecmp(argv[narg], "--use-prisms")) use_prisms = true;
for (int narg = 1; narg < argc - 1; narg++)
if (!strcasecmp(argv[narg], "--GDSIIFile")) GDSIIFile = argv[narg + 1];
if (GDSIIFile != 0 && use_prisms == true)
fprintf(stderr, "warning: --use-prisms is ignored if --GDSIIFile is specified\n");
bend_flux(true, GDSIIFile, use_prisms);
bend_flux(false, GDSIIFile, use_prisms);
// success if we made it here
return 0;
}