Dear All,
I was checking convergency of 2D bragg example from tutorial and found out
that adding "s.set_epsilon(eps)" to the code
make convergency better, due to different subpixel anisotropic averaging. Is it
possible to do such averaging by default in
"structure" constructor?
Results: Q-factor with set_epsilon and without vs N of grid points in unit.
N with without time,s
10 -668.325 -342.898 ~20
20 -699.382 -633.049 ~45
40 -695.826 -750.176 ~160
80 -694.357 -744.446 ~1 500
160 -693.971 -723.766 ~9 000
Full code, to plot table above I was changing lines with many #.
#include <meep.hpp>
using namespace meep;
const double half_cavity_width = 0.5*0.6, air_slit_width = 0.4,
grating_periodicity = 0.6,
half_waveguide_width = 1.0,
num_air_slits = 5.0,
high_dielectric = 4.0, low_dielectric = 2.0;
const double pml_thickness = 1.0;
const double x_center = 7.7 + pml_thickness;
const double y_center = 3.5 + pml_thickness;
double eps(const vec &rr) {
// Displacement from center of cavity is r:
const vec r = rr - vec(x_center, y_center);
// First the air slits:
double dx = fabs(r.x()) - half_cavity_width;
if (dx < num_air_slits*grating_periodicity && dx > 0.0) {
while (dx > grating_periodicity) dx -= grating_periodicity;
if (dx < air_slit_width) return 1.0;
}
// Now check if the y value is within the waveguide:
if (fabs(r.y()) < half_waveguide_width) return high_dielectric;
// Otherwise we must be in the surrounding low dielectric:
return low_dielectric;
}
int main(int argc, char *argv[]) {
initialize mpi(argc, argv);
deal_with_ctrl_c();
const double amicron = 10; // ############## a micron is this many grid
points.
grid_volume vol = vol2d(2*x_center, 2*y_center, amicron);
const symmetry S = mirror(Y, vol) + rotate2(Y, vol);
structure s(vol, eps, pml(pml_thickness), S);
s.set_epsilon(eps); // ############################# new subpixel
averaging
const char *dirname = make_output_directory(__FILE__);
s.set_output_directory(dirname);
fields f(&s);
const double wavelength = 1.72;
const double freq = 1.0/wavelength;
f.add_point_source(Hy, freq, 5.0, 0.0, 5.0, vec(x_center,y_center));
while (f.time() < f.last_source_time() && !interrupt) f.step();
f.output_hdf5(Dielectric, vol.surroundings());
while (f.time() < 400.0 && !interrupt) f.step();
double next_print_time = 500.0;
monitor_point *p = NULL;
FILE *myout = create_output_file(dirname, "hy");
while (f.time() <= 2000.0 && !interrupt) {
// Now we'll start taking data!
f.step();
p = f.get_new_point(vec(x_center,y_center), p);
master_fprintf(myout, "%g\t%g\t%g\n", f.time(),
real(p->get_component(Hy)),
imag(p->get_component(Hy)));
if (f.time() >= next_print_time) {
f.output_hdf5(Hy, vol.surroundings());
master_printf("Energy is %g at time %g\n",
f.total_energy(), f.time());
next_print_time += 500.0;
}
}
master_fclose(myout);
complex<double> *amp, *freqs;
int num;
FILE *myfreqs = create_output_file(dirname, "freqs");
//master_fprintf(myfreqs, "\\begin{verbatim}\n");
master_printf("Harminving Hy...\n");
interrupt = 0; // Harminv even if we were interrupted.
p->harminv(Hy, &, &freqs, &num, 0.8*freq, 1.2*freq, 5);
for (int i=0;i<num;i++) {
master_fprintf(myfreqs, "%g\t%g\t%g\t%g\t%g\n",
real(freqs[i]), imag(freqs[i]),
-real(freqs[i])/imag(freqs[i]),
real(amp[i]), imag(amp[i]));
}
// master_fprintf(myfreqs, "%cend{verbatim}\n", '\\');
master_fclose(myfreqs);
// delete[] dirname;
}
_______________________________________________
meep-discuss mailing list
[email protected]
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
