From 63bd6deb5b7752484cbd0f41f14512ece8ed4d32 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Wed, 16 May 2018 22:05:56 -0400 Subject: [PATCH] WIP: mode-solver improvements (#333) * mpb match-freq should only search k in direction d * adjust for bloch-periodic boundaries in get_eigenmode, so that it works for both sources and mode decomposition * allow get_eigenmode to return NULL if mode was not found, clean up a memory leak in mode coefficients * whoops * add parity argument to get_eigenmode_coefficients (fises #294) * add missing get_eigenmode_coefficients parameters, fix defaults --- .../Python_Tutorials/Mode_Decomposition.md | 22 +-- python/Makefile.am | 1 + python/simulation.py | 7 +- src/meep.hpp | 21 ++- src/mpb.cpp | 151 +++++++++--------- 5 files changed, 107 insertions(+), 95 deletions(-) diff --git a/doc/docs/Python_Tutorials/Mode_Decomposition.md b/doc/docs/Python_Tutorials/Mode_Decomposition.md index 44d7f2b0f..b1357cfaf 100644 --- a/doc/docs/Python_Tutorials/Mode_Decomposition.md +++ b/doc/docs/Python_Tutorials/Mode_Decomposition.md @@ -26,25 +26,25 @@ def main(args): Lt = args.Lt # taper length Si = mp.Medium(epsilon=12.0) - + dair = 3.0 dpml = 3.0 - + sx = dpml + Lw + Lt + Lw + dpml sy = dpml + dair + w2 + dair + dpml cell_size = mp.Vector3(sx, sy, 0) - geometry = [ mp.Block(material=Si, center=mp.Vector3(0,0,0), size=mp.Vector3(mp.inf,w1,mp.inf)), + geometry = [ mp.Block(material=Si, center=mp.Vector3(0,0,0), size=mp.Vector3(mp.inf,w1,mp.inf)), mp.Block(material=Si, center=mp.Vector3(0.5*sx-0.5*(Lt+Lw+dpml),0,0), size=mp.Vector3(Lt+Lw+dpml,w2,mp.inf)) ] # form linear taper hh = w2 ww = 2*Lt - + # taper angle (CCW, relative to +X axis) rot_theta = math.atan(0.5*(w2-w1)/Lt) - + pvec = mp.Vector3(-0.5*sx+dpml+Lw,0.5*w1,0) cvec = mp.Vector3(-0.5*sx+dpml+Lw+0.5*ww,0.5*hh+0.5*w1,0) rvec = cvec-pvec @@ -69,7 +69,7 @@ def main(args): # mode frequency fcen = 0.15 - + sources = [ mp.EigenModeSource(src=mp.GaussianSource(fcen, fwidth=0.5*fcen), component=mp.Ez, size=mp.Vector3(0,sy-2*dpml,0), @@ -78,9 +78,9 @@ def main(args): eig_parity=mp.ODD_Z, eig_kpoint=mp.Vector3(0.4,0,0), eig_resolution=32) ] - + symmetries = [ mp.Mirror(mp.Y,+1) ] - + sim = mp.Simulation(resolution=resolution, cell_size=cell_size, boundary_layers=boundary_layers, @@ -90,7 +90,7 @@ def main(args): xm = -0.5*sx + dpml + 0.5*Lw; # x-coordinate of monitor mflux = sim.add_eigenmode(fcen, 0, 1, mp.FluxRegion(center=mp.Vector3(xm,0), size=mp.Vector3(0,sy-2*dpml))); - + sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mp.Vector3(xm,0), 1e-10)) bands = np.array([1],dtype=np.int32) # indices of modes for which to compute expansion coefficients @@ -98,13 +98,13 @@ def main(args): alpha = np.zeros(2*num_bands,dtype=np.complex128) # preallocate array to store coefficients vgrp = np.zeros(num_bands,dtype=np.float64) # also store mode group velocities mvol = mp.volume(mp.vec(xm,-0.5*sy+dpml),mp.vec(xm,+0.5*sy-dpml)) - sim.fields.get_eigenmode_coefficients(mflux, mp.X, mvol, bands, alpha, vgrp) + sim.fields.get_eigenmode_coefficients(mflux, mp.X, mvol, bands, mp.ODD_Z, alpha, vgrp) alpha0Plus = alpha[2*0 + 0]; # coefficient of forward-traveling fundamental mode alpha0Minus = alpha[2*0 + 1]; # coefficient of backward-traveling fundamental mode print("refl:, {}, {:.8f}".format(Lt, abs(alpha0Minus)**2)) - + if __name__ == '__main__': parser = argparse.ArgumentParser() parser.add_argument('-Lt', type=float, default=3.0, help='taper length (default: 3.0)') diff --git a/python/Makefile.am b/python/Makefile.am index 96e49886c..71e6e8493 100644 --- a/python/Makefile.am +++ b/python/Makefile.am @@ -1,4 +1,5 @@ HPPFILES= \ + $(top_srcdir)/src/meep.hpp \ $(top_srcdir)/src/meep_internals.hpp \ $(top_srcdir)/src/bicgstab.hpp \ $(top_srcdir)/src/meep/vec.hpp \ diff --git a/python/simulation.py b/python/simulation.py index 87f52cea8..3eebea1df 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -1024,13 +1024,16 @@ def get_dft_array(self, dft_obj, component, num_freq): else: raise ValueError("Invalid type of dft object: {}".format(dft_obj)) - def get_eigenmode_coefficients(self, flux, bands, kpoint_func=None): + def get_eigenmode_coefficients(self, flux, bands, eig_parity=mp.NO_PARITY, + eig_vol=None, eig_resolution=0, eig_tolerance=1e-7, kpoint_func=None): if self.fields is None: raise ValueError("Fields must be initialized before calling get_eigenmode_coefficients") + if eig_vol is None: + eig_vol = flux.where num_bands = len(bands) coeffs = np.zeros(2 * num_bands, dtype=np.complex128) - self.fields.get_eigenmode_coefficients(flux, np.array(bands, dtype=np.intc), coeffs, None, kpoint_func) + self.fields.get_eigenmode_coefficients(flux, eig_vol, np.array(bands, dtype=np.intc), eig_parity, eig_resolution, eig_tolerance, coeffs, None, kpoint_func) return coeffs diff --git a/src/meep.hpp b/src/meep.hpp index 04e259ad9..8e7130705 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1441,10 +1441,10 @@ class fields { // values of field components at arbitrary points in space. // call destroy_eigenmode_data() to deallocate it when finished. void *get_eigenmode(double omega_src, direction d, const volume where, - const volume eig_vol, int band_num, - const vec &kpoint, bool match_frequency=true, - int parity=0, double resolution=0.0, - double eigensolver_tol=1.0e-7, bool verbose=false); + const volume eig_vol, int band_num, + const vec &kpoint, bool match_frequency, + int parity, double resolution, + double eigensolver_tol, bool verbose=false); void add_eigenmode_source(component c, const src_time &src, direction d, const volume &where, @@ -1457,11 +1457,22 @@ class fields { std::complex A(const vec &) = 0); void get_eigenmode_coefficients(dft_flux flux, - int *bands, int num_bands, + const volume &eig_vol, + int *bands, int num_bands, int parity, + double eig_resolution, double eigensolver_tol, std::complex *coeffs, double *vgrp, kpoint_func user_kpoint_func=0, void *user_kpoint_data=0); + void get_eigenmode_coefficients(dft_flux flux, + int *bands, int num_bands, int parity, + std::complex *coeffs, + double *vgrp, kpoint_func user_kpoint_func=0, + void *user_kpoint_data=0) { + get_eigenmode_coefficients(flux, flux.where, bands, num_bands, parity, + 0.0, 1e-7, coeffs, vgrp, user_kpoint_func, user_kpoint_data); + } + // initialize.cpp: void initialize_field(component, std::complex f(const vec &)); void initialize_with_nth_te(int n); diff --git a/src/mpb.cpp b/src/mpb.cpp index 231328b24..df7a74f21 100644 --- a/src/mpb.cpp +++ b/src/mpb.cpp @@ -206,7 +206,7 @@ void *fields::get_eigenmode(double omega_src, direction d, const volume where, const volume eig_vol, int band_num, - const vec &kpoint, bool match_frequency, + const vec &_kpoint, bool match_frequency, int parity, double resolution, double eigensolver_tol, @@ -216,15 +216,23 @@ void *fields::get_eigenmode(double omega_src, /*- part 1: preliminary setup for calling MPB -----------------*/ /*--------------------------------------------------------------*/ + // if the mode region extends over the full computational grid and we are bloch-periodic + // in any direction, set the corresponding component of the eigenmode initial-guess + // k-vector to be the (real part of the) bloch vector in that direction. + vec kpoint(_kpoint); + LOOP_OVER_DIRECTIONS(v.dim, d) + if (float(eig_vol.in_direction(d) == float(v.in_direction(d)))) + if (boundaries[High][d]==Periodic && boundaries[Low][d]==Periodic) + kpoint.set_direction(d, real(k[d])); + //bool verbose=true; if (resolution <= 0.0) resolution = 2 * gv.a; // default to twice resolution int n[3], local_N, N_start, alloc_N, mesh_size[3] = {1,1,1}; - mpb_real k[3] = {0,0,0}, kcart[3] = {0,0,0}; + mpb_real k[3] = {0,0,0}; double s[3] = {0,0,0}, o[3] = {0,0,0}; mpb_real R[3][3] = {{0,0,0},{0,0,0},{0,0,0}}; mpb_real G[3][3] = {{0,0,0},{0,0,0},{0,0,0}}; mpb_real kdir[3] = {0,0,0}; - double kscale = 1.0; double match_tol = eigensolver_tol * 10; if (d == NO_DIRECTION || coordinate_mismatch(gv.dim, d)) @@ -267,46 +275,39 @@ void *fields::get_eigenmode(double omega_src, if (!quiet && verbose) master_printf("KPOINT: %g, %g, %g\n", k[0], k[1], k[2]); - // if match_frequency is true, all we need is a direction for k - // and a crude guess for its value; we must supply this if k==0. - if (match_frequency && k[0] == 0 && k[1] == 0 && k[2] == 0) { - k[d-X] = omega_src * sqrt(get_eps(eig_vol.center())); - if(!quiet && verbose) - master_printf("NEW KPOINT: %g, %g, %g\n", k[0], k[1], k[2]); - if (s[d-X] > 0) { - k[d-X] *= s[d-X]; // put k in G basis (inverted when we compute kcart) - if (fabs(k[d-X]) > 0.4) // ensure k is well inside the Brillouin zone - k[d-X] = k[d-X] > 0 ? 0.4 : -0.4; - if(!quiet && verbose) - master_printf("NEWER KPOINT: %g, %g, %g\n", k[0], k[1], k[2]); - } - } - for (int i = 0; i < 3; ++i) { n[i] = int(resolution * s[i] + 0.5); if (n[i] == 0) n[i] = 1; R[i][i] = s[i] = s[i] == 0 ? 1 : s[i]; G[i][i] = 1 / R[i][i]; // recip. latt. vectors / 2 pi } - for (int i = 0; i < 3; ++i) - for (int j = 0; j < 3; ++j) - kcart[i] += G[j][i] * k[j]; - double klen0 = sqrt(k[0]*k[0]+k[1]*k[1]+k[2]*k[2]); - double klen = sqrt(kcart[0]*kcart[0]+kcart[1]*kcart[1]+kcart[2]*kcart[2]); - if (klen == 0.0) { - if (match_frequency) abort("need nonzero kpoint guess to match frequency"); - klen = 1; - } - kdir[0] = kcart[0] / klen; - kdir[1] = kcart[1] / klen; - kdir[2] = kcart[2] / klen; - maxwell_data *mdata = create_maxwell_data(n[0], n[1], n[2], &local_N, &N_start, &alloc_N, band_num, band_num); if (local_N != n[0] * n[1] * n[2]) abort("MPI version of MPB library not supported"); + meep_mpb_eps_data eps_data; + eps_data.s = s; eps_data.o = o; eps_data.dim = gv.dim; eps_data.f = this; + set_maxwell_dielectric(mdata, mesh_size, R, G, meep_mpb_eps,NULL, &eps_data); + if (check_maxwell_dielectric(mdata, 0)) + abort("invalid dielectric function for MPB"); + + double kmatch = G[d-X][d-X] * k[d-X]; // k[d] in cartesian + kdir[d-X] = 1; // kdir = unit vector in d direction + + // if match_frequency is true, we need at least a crude guess for kmatch; + // which we automatically pick if kmatch == 0. + if (match_frequency && kmatch == 0) { + vec cen = eig_vol.center(); + kmatch = omega_src * sqrt(get_eps(cen)*get_mu(cen)); + k[d-X] = kmatch * R[d-X][d-X]; // convert to reciprocal basis + if (eig_vol.in_direction(d) > 0 && fabs(k[d-X]) > 0.4) // ensure k is well inside the Brillouin zone + k[d-X] = k[d-X] > 0 ? 0.4 : -0.4; + if(!quiet && verbose) + master_printf("NEW KPOINT: %g, %g, %g\n", k[0], k[1], k[2]); + } + set_maxwell_data_parity(mdata, parity); update_maxwell_data_k(mdata, k, G[0], G[1], G[2]); @@ -317,12 +318,6 @@ void *fields::get_eigenmode(double omega_src, abort("zero-frequency bands at k=0 are ill-defined"); } - meep_mpb_eps_data eps_data; - eps_data.s = s; eps_data.o = o; eps_data.dim = gv.dim; eps_data.f = this; - set_maxwell_dielectric(mdata, mesh_size, R, G, meep_mpb_eps,NULL, &eps_data); - if (check_maxwell_dielectric(mdata, 0)) - abort("invalid dielectric function for MPB"); - evectmatrix H = create_evectmatrix(n[0] * n[1] * n[2], 2, band_num, local_N, N_start, alloc_N); for (int i = 0; i < H.n * H.p; ++i) { @@ -345,10 +340,12 @@ void *fields::get_eigenmode(double omega_src, maxwell_zero_k_constraint, (void *) mdata); - mpb_real knew[3]; for (int i = 0; i < 3; ++i) knew[i] = k[i]; - mpb_real vgrp; // Re( W[0]* (dTheta/dk) W[0] ) = group velocity + // track #times change in kmatch increases to detect non-convergence + double dkmatch_prev = kmatch; + int count_dkmatch_increase = 0; + /*--------------------------------------------------------------*/ /*- part 2: newton iteration loop with call to MPB on each step */ /*- until eigenmode converged to requested tolerance */ @@ -367,7 +364,7 @@ void *fields::get_eigenmode(double omega_src, (am_master() && verbose && !quiet ? EIGS_VERBOSE : 0)); if (!quiet) master_printf("MPB solved for omega_%d(%g,%g,%g) = %g after %d iters\n", - band_num, knew[0],knew[1],knew[2], + band_num, k[0],k[1],k[2], sqrt(eigvals[band_num-1]), num_iters); if (match_frequency) { @@ -375,9 +372,9 @@ void *fields::get_eigenmode(double omega_src, evectmatrix_resize(&W[0], 1, 0); evectmatrix_resize(&W[1], 1, 0); for (int i = 0; i < H.n; ++i) - W[0].data[i] = H.data[H.p-1 + i * H.p]; + W[0].data[i] = H.data[H.p-1 + i * H.p]; - // compute the group velocity in the k direction + // compute the group velocity in the kdir direction maxwell_ucross_op(W[0], W[1], mdata, kdir); // W[1] = (dTheta/dk) W[0] mpb_real vscratch; evectmatrix_XtY_diag_real(W[0], W[1], &vgrp, &vscratch); @@ -388,20 +385,22 @@ void *fields::get_eigenmode(double omega_src, evectmatrix_resize(&W[1], band_num, 0); // update k via Newton step - kscale = kscale - (sqrt(eigvals[band_num - 1]) - omega_src) / (vgrp*klen0); + double dkmatch = (sqrt(eigvals[band_num - 1]) - omega_src) / vgrp; + kmatch = kmatch - dkmatch; if (!quiet && verbose) - master_printf("Newton step: group velocity v=%g, kscale=%g\n", vgrp, kscale); - if (kscale < 0 || kscale > 100) - abort("Newton solver not converging -- need a better starting kpoint"); - for (int i = 0; i < 3; ++i) knew[i] = k[i] * kscale; - update_maxwell_data_k(mdata, knew, G[0], G[1], G[2]); + master_printf("Newton step: group velocity v=%g, kmatch=%g\n", vgrp, kmatch); + count_dkmatch_increase += fabs(dkmatch) > fabs(dkmatch_prev); + if (count_dkmatch_increase > 4) + return NULL; + k[d-X] = kmatch * R[d-X][d-X]; + update_maxwell_data_k(mdata, k, G[0], G[1], G[2]); } } while (match_frequency && fabs(sqrt(eigvals[band_num - 1]) - omega_src) > omega_src * match_tol); if (!match_frequency) - omega_src = sqrt(eigvals[band_num - 1]); + omega_src = sqrt(eigvals[band_num - 1]); // cleanup temporary storage delete[] eigvals; @@ -474,9 +473,9 @@ void *fields::get_eigenmode(double omega_src, edata->s[0] = s[0]; edata->s[1] = s[1]; edata->s[2] = s[2]; - edata->k[0] = knew[0]; - edata->k[1] = knew[1]; - edata->k[2] = knew[2]; + edata->k[0] = k[0]; + edata->k[1] = k[1]; + edata->k[2] = k[2]; edata->center = eig_vol.center() - where.center(); edata->amp_func = default_amp_func; edata->band_num = band_num; @@ -537,7 +536,6 @@ void fields::add_eigenmode_source(component c0, const src_time &src, double resolution, double eigensolver_tol, complex amp, complex A(const vec &)) { - /*--------------------------------------------------------------*/ /* step 1: call MPB to compute the eigenmode */ /*--------------------------------------------------------------*/ @@ -549,6 +547,9 @@ void fields::add_eigenmode_source(component c0, const src_time &src, parity, resolution, eigensolver_tol); + if (!global_eigenmode_data) + abort("eigenmode solver failed to find the requested mode; you may need to supply a better guess for k"); + global_eigenmode_data->amp_func = A ? A : default_amp_func; src_time *src_mpb = src.clone(); @@ -582,9 +583,6 @@ void fields::add_eigenmode_source(component c0, const src_time &src, destroy_eigenmode_data( (void *)global_eigenmode_data); } -bool equal_float(double d1, double d2) -{ return ((float)d1)==((float)d2); } - /***************************************************************/ /* get eigenmode coefficients for all frequencies in flux */ /* and all band indices in the caller-populated bands array. */ @@ -601,7 +599,9 @@ bool equal_float(double d1, double d2) /* bands[nb] is stored in vgrp[nb*num_freqs + nf]. */ /***************************************************************/ void fields::get_eigenmode_coefficients(dft_flux flux, - int *bands, int num_bands, + const volume &eig_vol, + int *bands, int num_bands, int parity, + double eig_resolution, double eigensolver_tol, std::complex *coeffs, double *vgrp, kpoint_func user_kpoint_func, void *user_kpoint_data) @@ -611,9 +611,6 @@ void fields::get_eigenmode_coefficients(dft_flux flux, int num_freqs = flux.Nfreq; direction d = flux.normal_direction; bool match_frequency = true; - int parity = 0; // NO_PARITY - double resolution = a; - double eig_tol = 1.0e-4; if (d==NO_DIRECTION) abort("cannot determine normal direction in get_eigenmode_coefficients"); @@ -621,15 +618,7 @@ void fields::get_eigenmode_coefficients(dft_flux flux, if (S.multiplicity() > 1) abort("symmetries are not yet supported in get_eigenmode_coefficients"); - // if the flux region extends over the full computational grid and we are bloch-periodic - // in any direction, set the corresponding component of the eigenmode initial-guess - // k-vector to be the (real part of the) bloch vector in that direction. otherwise just - // set the initial-guess k component to zero. - vec kpoint(0.0,0.0,0.0); - LOOP_OVER_DIRECTIONS(this->v.dim, d) - if ( equal_float(flux.where.in_direction(d), this->v.in_direction(d) ) ) - if (boundaries[High][d]==Periodic && boundaries[Low][d]==Periodic) - kpoint.set_direction(d, real(this->k[d])); + vec kpoint(0.0,0.0,0.0); // default guess // loop over all bands and all frequencies for(int nb=0; nb0.0 ? 0 : 1) ] = cplus/sqrt(abs(normfac)); coeffs[ 2*nb*num_freqs + 2*nf + (vg>0.0 ? 1 : 0) ] = cminus/sqrt(abs(normfac)); + destroy_eigenmode_data((void*) mode_data); } } @@ -668,10 +662,10 @@ void fields::get_eigenmode_coefficients(dft_flux flux, /**************************************************************/ #else // #ifdef HAVE_MPB void *fields::get_eigenmode(double omega_src, - direction d, const volume where, - const volume eig_vol, - int band_num, - const vec &kpoint, bool match_frequency, + direction d, const volume where, + const volume eig_vol, + int band_num, + const vec &kpoint, bool match_frequency, int parity, double resolution, double eigensolver_tol, bool verbose) { @@ -700,12 +694,15 @@ void fields::add_eigenmode_source(component c0, const src_time &src, } void fields::get_eigenmode_coefficients(dft_flux flux, - int *bands, int num_bands, + const volume &eig_vol, + int *bands, int num_bands, int parity, + double eig_resolution, double eigensolver_tol, std::complex *coeffs, double *vgrp, kpoint_func user_kpoint_func, void *user_kpoint_data) -{ (void) flux; (void) bands; (void)num_bands; +{ (void) flux; (void) eig_vol; (void) bands; (void)num_bands; + (void) parity; (void) eig_resolution; (void) eigensolver_tol; (void) coeffs; (void) vgrp; (void) user_kpoint_func; (void) user_kpoint_data; abort("Meep must be configured/compiled with MPB for get_eigenmode_coefficient"); }