Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

WIP: Allow unequal DFT frequencies #1141

Closed
wants to merge 37 commits into from

Conversation

danielwboyce
Copy link

@danielwboyce danielwboyce commented Mar 2, 2020

Fixes #1070. In progress of changing the relevant DFT classes' properties so that they have an array of frequencies, not freq_min or freq_max. The existing constructors will be updated so that DFT objects can be defined with freq_min, freq_max, and Nfreq, but new constructors will also be added that use only freqs and Nfreq.

@smartalecH
Copy link
Collaborator

Some swig advice:

When wrapping overloaded C constructors, Swig needs all of the parameters to be independently type-checked (see docs). This is good practice anyway... If you want to check a numpy array called freqs, you can do something like this:

%typecheck(SWIG_TYPECHECK_POINTER, fragment="NumPy_Fragments") double* freqs {
    $1 = is_array($input);
}

(all within meep.i of course)

You'll then have to define a custom typemap. You can do something like this:

%typemap(in, fragment="NumPy_Macros") double* freqs {
    $1 = (double *)array_data($input);
}

Make sure you do this for all the fields."flux" methods (i.e. add_dft_flux, etc.)

You'll have to go through simulation.py and update all of the flux, force, far_field, eigenmode, etc methods.

@danielwboyce
Copy link
Author

As of now, the source code seems to be properly compiling (none of the constructors seem to have been broken). Now I'll need to update python/meep.i and scheme/meep.i to have support in the scripting languages, which should also make testing easier.

@danielwboyce
Copy link
Author

As of now the package compiles all the way through, but I'm sure I'm getting SWIG wrong. The version of python/meep.i that's in place has the minimal changes in place required to build the package with the changes to the C code and simulation.py. I notice that for some reason the overloaded constructors don't need me to define new typemaps or typechecks to compile; is it likely that the automatically-generated wrappers are correct?

While attempting to run the python/examples/mie_scattering.py tutorial example I'm getting the following error

Traceback (most recent call last):
  File "/home/tests/mie_scattering.testing_v.pass_frq_min_frq_max.py", line 51, in <module>
    sim.run(until_after_sources=10)
  File "/usr/local/lib/python3.6/site-packages/meep/simulation.py", line 2259, in run
    self._evaluate_dft_objects()
  File "/usr/local/lib/python3.6/site-packages/meep/simulation.py", line 1627, in _evaluate_dft_objects
    dft.swigobj = dft.func(*dft.args)
  File "/usr/local/lib/python3.6/site-packages/meep/simulation.py", line 1851, in _add_flux
    return self._add_fluxish_stuff(self.fields.add_dft_flux, freqs, nfreq, fluxes)
  File "/usr/local/lib/python3.6/site-packages/meep/simulation.py", line 1987, in _add_fluxish_stuff
    stuff = add_dft_stuff(vol_list, freqs, nfreq, *args)
  File "/usr/local/lib/python3.6/site-packages/meep/__init__.py", line 3959, in add_dft_flux
    return _meep.fields_add_dft_flux(self, *args)
NotImplementedError: Wrong number or type of arguments for overloaded function 'fields_add_dft_flux'.
  Possible C/C++ prototypes are:
    meep::fields::add_dft_flux(meep::volume_list const *,double *,int,bool)
    meep::fields::add_dft_flux(meep::volume_list const *,double *,int)
    meep::fields::add_dft_flux(meep::volume_list const *,double,double,int,bool)
    meep::fields::add_dft_flux(meep::volume_list const *,double,double,int)
    meep::fields::add_dft_flux(meep::direction,meep::volume const &,double *,int,bool)
    meep::fields::add_dft_flux(meep::direction,meep::volume const &,double *,int)
    meep::fields::add_dft_flux(meep::direction,meep::volume const &,double,double,int,bool)
    meep::fields::add_dft_flux(meep::direction,meep::volume const &,double,double,int)

Which seems to tell me that the add_dft_flux function is not playing nice with SWIG. Thing is there's no add_dft_flux function mentioned in python/meep.i! Was the wrapper getting generated automatically before? Will I need to make a typemap suite for it now? Which parameters will I need to check? All of them, or just the ones that differ between the overloaded versions of the function (double freq_min, double freq_max, double *freqs)?

@smartalecH
Copy link
Collaborator

The python/swig UI regarding DFT stuff is a little convoluted...largely because of the load balancing requirements. But luckily I think you only need some minor changes to get your code working.

First, let's remember that all the DFT functions using fluxes will call add_dft_flux through _add_fluxish_stuff:

meep/python/simulation.py

Lines 1820 to 1823 in 6c10e74

def _add_flux(self, fcen, df, nfreq, fluxes):
if self.fields is None:
self.init_sim()
return self._add_fluxish_stuff(self.fields.add_dft_flux, fcen, df, nfreq, fluxes)

You can see that _add_fluxish_stuff assumes one particular constructor form for the corresponding add_dft_stuff function (add_dft_flux in this case):

meep/python/simulation.py

Lines 1939 to 1955 in 6c10e74

def _add_fluxish_stuff(self, add_dft_stuff, fcen, df, nfreq, stufflist, *args):
vol_list = None
for s in stufflist:
v = Volume(center=s.center, size=s.size, dims=self.dimensions,
is_cylindrical=self.is_cylindrical)
d0 = s.direction
d = self.fields.normal_direction(v.swigobj) if d0 < 0 else d0
c = mp.direction_component(mp.Sx, d)
v2 = Volume(center=s.center, size=s.size, dims=self.dimensions,
is_cylindrical=self.is_cylindrical).swigobj
vol_list = mp.make_volume_list(v2, c, s.weight, vol_list)
stuff = add_dft_stuff(vol_list, fcen - df / 2, fcen + df / 2, nfreq, *args)
vol_list.__swig_destroy__(vol_list)
return stuff

This is problematic for two reasons:

  1. First, you need to generalize the python code so we can either use the original constructor (i.e. specify center, width, and the number of freqs) or pass in an array (like you've already implemented on your branch).
  2. You need to explicitly typecheck and typemap every argument of all the possible constructors. Even if it compiles, that doesn't mean it works! The typechecking mechanism is how swig determines which constructor to use. It assigns a weight to each paramter based on it's type and maps the cumulative weight to a table. You can see more about that in the docs I posted above.

Feel free to check out PR #1095 for an example. I "overloaded" add_dft_fields (in reality I just added a default argument, but SWIG treats this as an overloaded constructor). I had to typecheck the remaining parameters. I was getting the same error as you until I realized this small nuance.

Looking forward to this feature!

@ChristopherHogan
Copy link
Contributor

ChristopherHogan commented Mar 12, 2020

I think what you want is to apply the predefined numpy typemaps. If you have double *freqs, int Nfreqs pairs of arguments for the new functions and they're all named consistently, all you need in python/meep.i is

%apply (double *IN_ARRAY1, int DIM1) {(double *freqs, int Nfreqs)};

@ChristopherHogan
Copy link
Contributor

self.dft_objects.append(flux)
return flux

def _add_flux(self, fcen, df, nfreq, fluxes):
def add_flux(self, fcen, df, nfreq, *fluxes):
freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq).tolist()
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You can remove the tolist() call since the typemaps are expecting a numpy array.

python/meep.i Outdated
$1 = is_array($input);
}

%typemap(in, fragment="NumPy_Macros") double* freqs {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This typemap isn't needed. It's covered by the %apply below.

python/meep.i Outdated
$1 = PyInteger_Check($input);
}

%typemap(in) int Nfreqs {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This typemap isn't needed. It's covered by the %apply below.

@danielwboyce
Copy link
Author

Removing tolist where extraneous and the extra typechecks doesn't seem to have fixed the problem.

src/dft_ldos.cpp Outdated
Nomega = Nfreq;
double om[Nomega];
for (int i = 0; i < Nfreq; i++) { om[i] = freqs[i] * 2 * pi; }
omegas = om;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This won't work. om will no longer exist once this scope exits. You need to allocate on the heap with e.g. new

src/dft_ldos.cpp Outdated
domega = 0;
Nomega = 1;
// TODO: overload dft_ldos constructor
dft_ldos::dft_ldos(double *freqs, int Nfreq) {
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If the object makes a copy of freqs, then this should be declared as const double *freqs.

src/dft.cpp Outdated
bool use_centered_grid, int vc) {
if (Nfreq < 1) abort("Nfreq must be at least 1");
double dfreq = (freq_max - freq_min) / (Nfreq - 1);
double freqs[Nfreq];
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same problem: this does not persist.

src/meep.hpp Outdated
int Nomega;
double *omegas;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

might want to use std::vector<double> here.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

also below

@stevengj
Copy link
Collaborator

Before you worry about getting the Python interface to work, I would to get the C++ tests to work, e.g. tests/bend-flux-ll.cpp and tests/dft-fields.cpp.

@stevengj
Copy link
Collaborator

In fact, let's separate the Python stuff into a separate PR. Just try to get the C++ plumbing in place in this PR.

@danielwboyce
Copy link
Author

Let's also add a C++ test for the new plumbing in this PR.

@oskooi
Copy link
Collaborator

oskooi commented Mar 23, 2020

Note that this PR is no longer necessary as it has been replaced by #1154.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

allow unequal DFT frequencies
5 participants