flatironinstitute / FMM3D

Flatiron Institute Fast Multipole Libraries --- This codebase is a set of libraries to compute N-body interactions governed by the Laplace and Helmholtz equations, to a specified precision, in three dimensions, on a multi-core shared-memory machine.
https://fmm3d.readthedocs.io
Other
91 stars 36 forks source link

Helmholtz FMM inaccurate when Re(zk)<0 #16

Closed bobbielf2 closed 2 years ago

bobbielf2 commented 3 years ago

Hi folks,

I was in MATLAB testing hfmm3d with a negative wavenumber zk<0 (see code below), but keep getting only ~2 digits when using a moderate tolerance (e.g. 1e-6). (The original goal was to do fast adjoint mat-vec.)

Here is some more info:

  1. Tests would fail when Re(zk)<0
  2. However, can get full 15 digits once the number of sources ns is small enough (e.g. ns=1000)
  3. Can also get 15 digits once tolerance tol is small enough (e.g. tol=1e-12)

So I am getting either 2 digits or 15 digits. Could you confirm if this is a bug or am I doing something wrong here?

Thanks, Bobbie

ns = 4000;
srcinfo.sources = rand(3,ns);
srcinfo.charges = rand(1,ns)+1i*rand(1,ns);
stmp = srcinfo.sources;
pg = 1;

zk = complex(-1.1);
tol = 1e-6;

fprintf('Re(zk) = %.1f\ntol = %.1e\n',real(zk),tol)

U1 = hfmm3d(tol,zk,srcinfo,pg);
U2 = h3ddir(zk,srcinfo,stmp,pg);
U2.pot = U2.pottarg;
err = norm(U1.pot-U2.pot)/norm(U2.pot);
if err<tol, fprintf('Test successful, '), 
else, fprintf('Test failed, '); end
fprintf('err = %.2e\n',err);
mrachh commented 3 years ago

Hi, Thanks for bringing this to our notice. Currently, the Helmholtz wave number zk should be in the first quadrant. In case you do want to run the calculation for zk in the second quadrant, you can conjugate the density on input and the result on output. We will update our documentation soon to reflect this restriction.

Regards, Manas

bobbielf2 commented 3 years ago

Thanks, Manas! I will use the double conjugation you said.

On Wed, Nov 18, 2020, 8:05 AM mrachh notifications@github.com wrote:

Hi, Thanks for bringing this to our notice. Currently, the Helmholtz wave number zk should be in the first quadrant. In case you do want to run the calculation for zk in the second quadrant, you can conjugate the density on input and the result on output. We will update our documentation soon to reflect this restriction.

Regards, Manas

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/flatironinstitute/FMM3D/issues/16#issuecomment-729699118, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEHNZ7BD6BL63KP47AGIVJLSQPIDHANCNFSM4TZQNCAQ .

ahbarnett commented 3 years ago

Hello - just wanted to chime in here that 2nd or 3rd quadrant are negative Re part of k, which is modified Helmholtz (decaying) rather than Helmholtz (oscillating). I'm sure you can't access that via conjugation (which would only get you to 4th quadrant). A different FMM is needed. (Unless your kernel scale is << system scale, when no FMM is needed, just a NN search.) Best, Alex

On Wed, Nov 18, 2020 at 10:27 AM bobbielf2 notifications@github.com wrote:

Thanks, Manas! I will use the double conjugation you said.

On Wed, Nov 18, 2020, 8:05 AM mrachh notifications@github.com wrote:

Hi, Thanks for bringing this to our notice. Currently, the Helmholtz wave number zk should be in the first quadrant. In case you do want to run the calculation for zk in the second quadrant, you can conjugate the density on input and the result on output. We will update our documentation soon to reflect this restriction.

Regards, Manas

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub < https://github.com/flatironinstitute/FMM3D/issues/16#issuecomment-729699118 , or unsubscribe < https://github.com/notifications/unsubscribe-auth/AEHNZ7BD6BL63KP47AGIVJLSQPIDHANCNFSM4TZQNCAQ

.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/flatironinstitute/FMM3D/issues/16#issuecomment-729754303, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACNZRSU5AHEN7GHSMUDSR63SQPRVZANCNFSM4TZQNCAQ .

-- *---------------------------------------------------------------------~^`^~._.~' |\ Alex H. Barnett Center for Computational Mathematics, Flatiron Institute | \ http://users.flatironinstitute.org/~ahb 646-876-5942

lgreengard commented 3 years ago

Hi - negative real part is actually the other side of the Fourier transform and annoying to incorporate in the code since it implies a different propagating direction, a different radiation condition, etc. Double conjugation gets around that. (Modified Helmholtz corresponds to pure imaginary k.)

On 11/18/20 2:58 PM, Alex Barnett wrote:

Hello - just wanted to chime in here that 2nd or 3rd quadrant are negative Re part of k, which is modified Helmholtz (decaying) rather than Helmholtz (oscillating). I'm sure you can't access that via conjugation (which would only get you to 4th quadrant). A different FMM is needed. (Unless your kernel scale is << system scale, when no FMM is needed, just a NN search.) Best, Alex

On Wed, Nov 18, 2020 at 10:27 AM bobbielf2 notifications@github.com wrote:

Thanks, Manas! I will use the double conjugation you said.

On Wed, Nov 18, 2020, 8:05 AM mrachh notifications@github.com wrote:

Hi, Thanks for bringing this to our notice. Currently, the Helmholtz wave number zk should be in the first quadrant. In case you do want to run the calculation for zk in the second quadrant, you can conjugate the density on input and the result on output. We will update our documentation soon to reflect this restriction.

Regards, Manas

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <

https://github.com/flatironinstitute/FMM3D/issues/16#issuecomment-729699118

, or unsubscribe <

https://github.com/notifications/unsubscribe-auth/AEHNZ7BD6BL63KP47AGIVJLSQPIDHANCNFSM4TZQNCAQ

.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub

https://github.com/flatironinstitute/FMM3D/issues/16#issuecomment-729754303, or unsubscribe

https://github.com/notifications/unsubscribe-auth/ACNZRSU5AHEN7GHSMUDSR63SQPRVZANCNFSM4TZQNCAQ .

-- *---------------------------------------------------------------------~^`^~._.~' |\ Alex H. Barnett Center for Computational Mathematics, Flatiron Institute | \ http://users.flatironinstitute.org/~ahb 646-876-5942

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/flatironinstitute/FMM3D/issues/16#issuecomment-729919175, or unsubscribe https://github.com/notifications/unsubscribe-auth/AANMHYPSSHBTUIA6O6OI6X3SQQROPANCNFSM4TZQNCAQ.

bobbielf2 commented 3 years ago

Howdy Alex! I am not familiar with modified Helmholtz. Are you talking about the Re(k^2)<0 case (i.e. screened Coulomb/Yukawa potential)? That would indeed require different FMM. To be more specific, what I was trying to do (conjugate transpose mat-vec) is still the unmodified Helmholtz. E.g. SLPmatrix'*x = conj(SLPmatrix*conj(x)) (the case with quadrature weights is analogous). Thanks for joining in (and thanks for your input on the other thread)! Best, Bobbie

bobbielf2 commented 3 years ago

Hi Leslie, Thanks for the explanation. (I see that hfmm3d works for pure imaginary k.) Best, Bobbie

ahbarnett commented 3 years ago

My bad - was thinking k^2 (Helmholtz parameter), not k. Best, Alex

On Wed, Nov 18, 2020 at 3:44 PM Leslie Greengard notifications@github.com wrote:

Hi - negative real part is actually the other side of the Fourier transform and annoying to incorporate in the code since it implies a different propagating direction, a different radiation condition, etc. Double conjugation gets around that. (Modified Helmholtz corresponds to pure imaginary k.)

  • Leslie

On 11/18/20 2:58 PM, Alex Barnett wrote:

Hello - just wanted to chime in here that 2nd or 3rd quadrant are negative Re part of k, which is modified Helmholtz (decaying) rather than Helmholtz (oscillating). I'm sure you can't access that via conjugation (which would only get you to 4th quadrant). A different FMM is needed. (Unless your kernel scale is << system scale, when no FMM is needed, just a NN search.) Best, Alex

On Wed, Nov 18, 2020 at 10:27 AM bobbielf2 notifications@github.com wrote:

Thanks, Manas! I will use the double conjugation you said.

On Wed, Nov 18, 2020, 8:05 AM mrachh notifications@github.com wrote:

Hi, Thanks for bringing this to our notice. Currently, the Helmholtz wave number zk should be in the first quadrant. In case you do want to run the calculation for zk in the second quadrant, you can conjugate the density on input and the result on output. We will update our documentation soon to reflect this restriction.

Regards, Manas

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <

https://github.com/flatironinstitute/FMM3D/issues/16#issuecomment-729699118

, or unsubscribe <

https://github.com/notifications/unsubscribe-auth/AEHNZ7BD6BL63KP47AGIVJLSQPIDHANCNFSM4TZQNCAQ

.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub

< https://github.com/flatironinstitute/FMM3D/issues/16#issuecomment-729754303 , or unsubscribe

< https://github.com/notifications/unsubscribe-auth/ACNZRSU5AHEN7GHSMUDSR63SQPRVZANCNFSM4TZQNCAQ

.

--

*---------------------------------------------------------------------~^`^~._.~' |\ Alex H. Barnett Center for Computational Mathematics, Flatiron Institute | \ http://users.flatironinstitute.org/~ahb 646-876-5942

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub < https://github.com/flatironinstitute/FMM3D/issues/16#issuecomment-729919175>,

or unsubscribe < https://github.com/notifications/unsubscribe-auth/AANMHYPSSHBTUIA6O6OI6X3SQQROPANCNFSM4TZQNCAQ .

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/flatironinstitute/FMM3D/issues/16#issuecomment-729942728, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACNZRSS2WGKIRZOPAJCMOVLSQQW4HANCNFSM4TZQNCAQ .

-- *---------------------------------------------------------------------~^`^~._.~' |\ Alex H. Barnett Center for Computational Mathematics, Flatiron Institute | \ http://users.flatironinstitute.org/~ahb 646-876-5942