Dear Konstantin,
if you call the structure constructor with these arguments it creates
still the structure, but without default subpixel averaging. However
you can set additional arguments to the constructor. Here is the
headline of the structure constructor, defined in meep.hpp (as I
remember...).
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);
In this case the solution for you would be to use the constructor like this:
structure s(vol, eps, pml(pml_thickness), S, 0, 0.5, true);
In contrast the set_epsilon function uses anisotropic subpixel
averaging by default.
I hope this helps,
Ors
2011/3/3 Konstantin Ladutenko <[email protected]>:
> 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
>
_______________________________________________
meep-discuss mailing list
[email protected]
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
