I've been programming the examples that are listed in libctl in C++ and on the third example I am getting different results with the flux output. I did specify the C++ program to use only real fields. The libctl code is on the meep site and the C++ I'll post, but it is a bit long.
Brian
/* -*- Mode: C; indent-tabs-mode: t; c-basic-offset: 4; tab-width: 4 -*- */ /* * main.cc * Copyright (C) Brian Beardall 2010 <[email protected]> * * meepsample2 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 3 of the License, or * (at your option) any later version. * * meepsample2 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, see <http://www.gnu.org/licenses/>. */ #include <iostream> #include <meep.hpp> #include <cstring> #include <cstdlib> using namespace std; using namespace meep; //#define SIZEX "sizex"; //#define SIZEY sizey; //#define PAD pad; //#define WIDTH width; #define TRUE 1 #define FALSE 0 const char *SIZEX = "sizex"; const char *SIZEY = "sizey"; const char *PAD = "pad"; const char *WIDTH = "width"; const char *BEND = "bend"; const char *NFREQ = "nfreq"; int sx = 0; int sy = 0; float pad = 0; float w = 0; int nfreq = 0; double eps_bend(const vec &p) { if (p.y() > pad && p.y() < pad+w && p.x() < sx-pad) return 12.0; if (p.y() > pad && p.x() > sx-pad-w && p.x() < sx-pad) return 12.0; return 1.0; } double eps_straight(const vec &p) { if (p.y() > pad && p.y() < pad+w) return 12.0; return 1.0; } int main(int argc, char **argv) { initialize mpi(argc, argv); // do this even for non-MPI Meep char *str, *argi; bool bend = FALSE; if (argc > 1) { for(int size_argc = 1; size_argc < argc; size_argc++) { int n,m,dashcount=0; for(n = 0;*(*(argv+size_argc)+n) != '\0' && *(*(argv+size_argc)+n) != '='; n++) { if(*(*(argv+size_argc)+n) == '-') { dashcount++; } } // cout << dashcount << " " << n << endl; str = new char [n+1]; if(dashcount == 2) { for (m = 0; m+dashcount < n; m++) { *(str+m) = *(*(argv+size_argc)+m+dashcount); } *(str+n) = '\0'; if(*(*(argv+size_argc)+n) == '=') { int arg_length; n +=1; for(arg_length=0;*(*(argv+size_argc)+arg_length+n) != '\0';arg_length++) { } if(arg_length > 0) { m = 0; argi = new char [arg_length]; for(m = 0;m+n < arg_length+n; m++) { *(argi+m) = *(*(argv+size_argc)+m+n); } } } // cout << argi << endl; // cout << "separate\n"; // cout << str << endl; if (!strcmp(str,SIZEX)) { sx = atoi(argi); delete argi; } else if(!strcmp(str,SIZEY)) { sy = atoi(argi); delete argi; } else if (!strcmp(str,PAD)) { pad = atof(argi); delete argi; } else if(!strcmp(str,WIDTH)) { w = atof(argi); delete argi; } else if(!strcmp(str,NFREQ)) { nfreq = atoi(argi); delete argi; } else if(!strcmp(str,BEND)) { bend = TRUE; } } delete str; } } cout << sx << " " << sy << " " << w << " " << pad << endl; double resolution = 10; // pixels per distance grid_volume v = vol2d(sx,sy, resolution); // 5x10 2d cell fields *f; structure *s; if(bend) { s = new structure (v, eps_bend, pml(1.0)); } else { s = new structure (v, eps_straight, pml(1.0)); } f = new fields (s); (*f).output_hdf5(Dielectric, v.surroundings()); //double freq = 1/(2*sqrt(2)); double freq = 0.15; double fwidth = 0.1; double fmin = freq-fwidth/2; double fmax = freq+fwidth/2; h5file *file = bend ? (*f).open_h5file("ez-bend",h5file::WRITE) : (*f).open_h5file("ez-straight",h5file::WRITE); h5file *fflux = bend ? (*f).open_h5file("refl-flux",h5file::READONLY) : (*f).open_h5file("refl-flux",h5file::WRITE); //continuous_src_time src(freq, 20,0,200); gaussian_src_time src(freq,fwidth); (*f).use_real_fields(); bend ? (*f).add_point_source(Ez, src, vec(1, w/2+pad)):(*f).add_point_source(Ez, src,vec(1,w/2+pad)); double twrite = (*f).time(); //cout << f.time() << endl; vec rt(1.5,pad+w/2-w),rb(1.5,pad+w/2+w); vec *it,*ib; if(bend) { it = new vec(sx-pad-w/2-w,sy-1.5), ib = new vec(sx-pad-w/2+w,sy-1.5); } else { it = new vec(sx-1.5,pad+w/2-w); ib = new vec(sx-1.5,pad+w/2+w); } vec *at; at = new vec(sx-pad-w/2-w,sy-1.5); vec ab(sx-1.5,pad+w/2+w); //it = new vec(sx-1.5,pad+w/2-w); ib = new vec(sx-1.5,pad+w/2+w); cout << (*it).x() << ", " << (*it).y() << ", " << (*it).z() << endl; cout << (*ib).x() << ", " << (*ib).y() << ", " << (*ib).z() << endl; cout << (*at).x() << ", " << (*at).y() << ", " << (*at).z() << endl; cout << ab.x() << ", " << ab.y() << ", " << ab.z() << endl; //flux_vol *right = (*f).add_flux_plane(*it,*ib); //flux_vol *left = (*f).add_flux_plane(rt,rb); volume plane(*it,*ib); volume refl(rt,rb); cout << "center: " << plane.center().x() << ", " << plane.center().y() << ", " << plane.center().z() << endl; // plane.dim = D3; cout << "diameter: " << plane.diameter() << ", " << plane.dim << endl; cout << "fields: " << (*f).v.center().x() << ", " << (*f).v.center().y() << ", " << (*f).v.center().z() << endl; cout << rt.x() << ", " << rt.y() << ", " << rt.z() << endl; cout << rb.x() << ", " << rb.y() << ", " << rb.z() << endl; direction d = (*f).normal_direction(plane); int c = direction_component(Sx,d); cout << d << endl; volume_list *vl = new volume_list(plane,c); dft_flux flux1 = (*f).add_dft_flux(vl,fmin,fmax,nfreq); cout << "added flux\n"; dft_flux flux2 = (*f).add_dft_flux_plane(refl,fmin,fmax,nfreq); if(bend) { flux2.load_hdf5(fflux); flux2.scale_dfts(-1.0); } while ((*f).time() < (*f).last_source_time()) { // f.output_hdf5(Ex, v.surroundings(),&file,1); if(twrite + 0.6 <= (*f).time()) { (*f).output_hdf5(Ez, v.surroundings(),file,1); //cout << flux1.flux() << endl; //cout << f.time() << " " << twrite+0.6 << endl; twrite = (*f).time(); } (*f).step(); } if(bend) { while (1e-3 < abs((*f).get_field(Ez,vec(sx - 1.5,pad+w/2)))) { // f.output_hdf5(Ex, v.surroundings(),&file,1); if(twrite + 0.6 <= (*f).time()) { (*f).output_hdf5(Ez, v.surroundings(),file,1); //cout << flux1.flux() << endl; //cout << f.time() << " " << twrite+0.6 << endl; twrite = (*f).time(); } (*f).step(); } } else { while (1e-3 < abs((*f).get_field(Ez,vec(sx-1.5,pad+w/2)))) { // f.output_hdf5(Ex, v.surroundings(),&file,1); if(twrite + 0.6 <= (*f).time()) { (*f).output_hdf5(Ez, v.surroundings(),file,1); //cout << flux1.flux() << endl; //cout << f.time() << " " << twrite+0.6 << endl; twrite = (*f).time(); } (*f).step(); } } double sfreq = freq - fwidth/2; double dfreq = fwidth/(nfreq-1); for (int i = 0; i < nfreq;i++) { //cout << "flux: " << sfreq << ", " << *(flux1.flux()+i) << ", " << *(flux2.flux()+i) << endl; master_printf("flux: %g, %g, %g\n",sfreq,*(flux1.flux()+i),*(flux2.flux()+i)); sfreq+=dfreq; } if(!bend) flux2.save_hdf5(fflux); //f.output_hdf5(Ez, v.surroundings()); // f.output_hdf5( delete it; delete ib; delete vl; delete s; delete f; return 0; }
_______________________________________________ meep-discuss mailing list [email protected] http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
