NanoComp / meep

free finite-difference time-domain (FDTD) software for electromagnetic simulations
GNU General Public License v2.0
1.2k stars 613 forks source link

Correct the weights for dipole sources in cylindrical coordinates #2476

Open smartalecH opened 1 year ago

smartalecH commented 1 year ago

As discussed, it's challenging to get the weights just right for a dipole source in cylindrical coordinates. This PR attempts to address that.

First, I modified sources.cpp to include the r factor from the dV0 and dV1 weight factors, just like what's implemented when processing DFT fields. Previously, this factor was completely omitted.

Second, for dipole sources, I include the 1/r term that is necessary for a proper delta function in cylindrical coordinates.

Finally, I added some convenience functions that sum up all the source amplitudes (i.e. perform a sort of "integration") to check if what we're doing is right.

Using these three changes, we should be able to confirm that for a dipole source, it's integral should always be 1 (modulo the resolution scale factor) even for dipoles that lie between grid points. Previously, that wasn't the case. Now, it is. I ran a simple experiment that sweeps a dipole between two grid points and sum up the corresponding weights. Indeed, it sums to 1 for all values of $r$ (there's some ambiguity here at $r=0$).

weights.

Currently, some of the tests are failing. This could be due to a few different factors. Before we tackle that, I want to ensure we get the weight formulation down first.

stevengj commented 1 year ago

It might be worthwhile to think this through more carefully from a finite-element perspective (thinking of finite differences as equivalent to finite elements with first-order elements and a uniform grid, which is typically the case). Then the normalization is an integral of the currents "tent functions" times r.

This also gives you a better prescription for what to do with an r=0 gridpoint, which in FEM would be "half a tent" function.

oskooi commented 1 year ago

It seems this PR does not produce the expected result for the test in #2108. The radiated flux normalized by the dipole position $r^2$ is not a constant value independent or $r$. Results for resolution of 50.

freespace_cyl_flux_vs_rpos_PR2476

The results do not change with increasing resolution. Results below for resolution of 100. This is similar to the master branch as demonstrated in #2470 (comment).

freespace_cyl_flux_vs_rpos_res100_PR2476