need printf statement to pass value in eels branch r1301

[GoogleCodeExporter]

file GenerateB.c

lines 465-466:
   465   ctemp = imExp(t3*z/t1)* t4;
   466   printf("\n");

Without line 466 ctemp value is a NAN (or just not passed? Cannot check, 
because if I print out the value, the statement is present, and the bug is not 
there).

Since it's calculated for each dipole, many lines get printed.

Original issue reported on code.google.com by t.ostase...@gmail.com on 22 Jan 2014 at 5:59

[GoogleCodeExporter]

I couldn't reproduce a bug. I've checked out the current branch state and 
compiled it as is and with printf commented out (using 'make seq 
OPTIONS=DEBUG'). Then I ran both with 
adda.exe -grid 4 -beam electron 10 0 0 0
and get the same answer (in particular, gamma=1.10992035e-013).

Then I compiled the code with 'make seq' and ran with the same command - again 
everything working fine.

The question is how do you catch (identify) the NAN value?

Also a side note - imExp expects real argument, so putting creal inside imExp 
is logical. This should be covered by default typecast (from complex to real), 
but explicit cast is better for maintainability. Alternatively (recommended), 
you may declare real variables as double (use some other names than t#) rather 
than complex.

Original comment by yurkin on 23 Jan 2014 at 5:36

[GoogleCodeExporter]

 > The question is how do you catch (identify) the NAN value?

The iterative solver shows residue as a NAN and stops calculations after the 
first iteration, returning overall calculated value as a NAN as well.

 > I couldn't reproduce a bug.

I could reproduce the bug on several different linux machines (including a big 
cluster (Darwin) that I assume is correctly configured and up to date), but 
haven't tried to do that on windows machines (I assume you are using one, 
looking at 'adda.exe'). I will try to make sense of that at some later point. 
Maybe changing compilers or their versions?

Also, the bug seems to only appear when the beam is close to the particle, 
which brings me to another question - I wasn't able to find any descriptions 
where the (automatically generated) particles are located (and if centered) in 
the lab frame and if it is consistent with changing particles. The information 
is very important, as the EELS measurements are very sensitive to particle - 
beam distance, hence setting beam location becomes difficult when you are not 
sure where for example, sphere begins and ends.

 > Alternatively (recommended), you may declare real variables as double (use some other names than t#) rather than complex

Thanks for clearing that up, I changed it this way, just not committed yet.

Original comment by t.ostase...@gmail.com on 29 Jan 2014 at 10:00

[GoogleCodeExporter]

I just checked, and apparently the bug is present when both:
- The beam is not R-symmetric ( checked by simply removing the beam symmetry 
changing lines in the initialization code, such that the beam is always 
maximally symmetric (first time done to temporarily prevent double calculations 
that always result in the same answer) )

- The beam is sufficiently far away (i.e -gird 4 -beam electron 300 4 0 0 runs 
fine, but -grid 4 -beam electron 300 3 0 0 returns this:

 > CoupleConstant:0.005259037197+1.843854148e-05i
 > x_0 = 0
 > RE_000 = -NAN
 > gamma   = -nan

Original comment by t.ostase...@gmail.com on 29 Jan 2014 at 10:13

[GoogleCodeExporter]

OK, I seem to understand what is going on. I have reproduced the bug with the 
command line:
adda.exe -grid 16 -size 8 -no_vol_cor -beam electron 300 3.25 0.25 0 -store_beam

The problem is that beam passes exactly through one or several dipoles, leading 
to infinite values of incident beam at these dipoles. This can be checked by 
looking at the produced files IncBeam-X,Y.

So the fundamental question is whether the EELS-DDA formulation is valid at 
all, when the beam intersects the particle. It is definitely breaking down for 
exact hit into the dipole, but I guess that it may also be incorrect for beam 
passing between the dipoles (due to some incorrect handling of integral over 
the singularity, etc.). For instance, for point-dipole incident field we 
decided that dipole should be located strictly outside of particle, and this is 
a user's responsibility to ensure that. But if you feel that such experiment is 
OK (and some other methods exist to simulate it), we need to think how to 
isolate the corresponding singularity in the DDA.

Concerning the reference frames (RFs) - see Section 6.1 "Reference frames" of 
the manual. In short, all RFs are centered at the center of the dipole 
computational grid. The command line specifies beam center and propagation 
direction in the laboratory reference frame, which are then automatically 
transformed into the particle RF (potentially, rotated with respect to the 
laboratory one). You can get the better feeling of particle RF by looking at 
the coordinates of the dipoles (first three columns) in, e.g., files IncBeam-#.

Original comment by yurkin on 30 Jan 2014 at 2:32

[GoogleCodeExporter]

It's been a while since I looked into it, as I have been otherwise busy, but I 
also get the nan value when the beam does not pass directly through a dipole:

./adda -grid 100 -size 0.08 -beam electron 200 0.09 0.0 0 -store_beam 
-no_vol_cor

still results in nan values for fields, although beam passes at (0.09, 0, 0) 
and dipoles are distributed in region (-0.04:0.04, -0.04:0.04, -0.04:0.04)

I will update the issue if I find anything related to it / fix it when I have 
the time

Thanks for the suggestion though

Original comment by t.ostase...@gmail.com on 14 Feb 2014 at 10:59

[GoogleCodeExporter]

That's strange indeed. I've just tried the same command line and everything 
went fine...

That is rather improbable, but the problem may be connected with compiler. 
Which compiler are you using? When you have time, please post first lines of 
the output of 'adda -V'. Also try recompiling adda with 'OPTIONS=DEBUGFULL'.

Original comment by yurkin on 14 Feb 2014 at 11:56

[GoogleCodeExporter]

Sorry for the delay, these are the relevant lines:

./adda -V
ADDA v.1.3b4 (r1329M)
Sequential version
Built with GNU compilers version 4.8.2 (64-bit)
Extra build options: DEBUGFULL

And the output of the same line is as follows:

./adda -grid 100 -size 0.08 -beam electron 200 0.09 0.0 0 -store_beam 
-no_vol_cor
DEBUG: ../ADDAmain.c:81: Reading command line finished 
DEBUG: ../ADDAmain.c:85: Initialization of shape finished 
all data is saved in 'run104_sphere_g100_m1.5'
DEBUG: ../comm.c:512: 0 :  0 100 100 1000000 200 

DEBUG: ../ADDAmain.c:94: Make particle finished 
box dimensions: 100x100x100
lambda: 6.283185307   Dipoles/lambda: 7853.98
Required relative residual norm: 1e-05
Total number of occupied dipoles: 523984
DEBUG: ../ADDAmain.c:98: Calculator started 
DEBUG: ../calculator.c:669: InitDmatrix started 
Memory usage for MatVec matrices: 282.0 MB
DEBUG: ../fft.c:1135: Initialize FFT (1st part) 
DEBUG: ../fft.c:581: FFTW library version: fftw-3.3.3-sse2
     compiler: gcc -std=gnu99 -march=x86-64 -mtune=generic -O2 -pipe -fstack-protector --param=ssp-buffer-size=4 -O3 -fomit-frame-pointer -malign-double -fstrict-aliasing -ffast-math
     codelet optimizations:  
Calculating Green's function (Dmatrix)
Fourier transform of Dmatrix......
Initializing FFTW3
DEBUG: ../calculator.c:671: InitDmatrix finished 
Total memory usage: 489.4 MB

here we go, calc Y

CoupleConstant:3.595029548e-11+8.616158303e-22i
DEBUG: ../CalculateE.c:876: Generating B 
Incident beam saved to file
DEBUG: ../CalculateE.c:881: Iterative solver started 
x_0 = E_inc
RE_000 = -NAN
DEBUG: ../CalculateE.c:883: Iterative solver finished 
DEBUG: ../CalculateE.c:497: Accumulating Eplane started 
DEBUG: ../CalculateE.c:500: Accumulating Eplane finished 
DEBUG: ../CalculateE.c:706: Calculation of EELS probability started 
gamma   = -nan
DEBUG: ../CalculateE.c:715: Calculation of EELS probability finished 

Original comment by t.ostase...@gmail.com on 24 Feb 2014 at 2:12

[GoogleCodeExporter]

When the incident beam is saved, for exactly the same configurations, lines are 
as follows (with the print("\n") line and without):

------------with-----------------
0.003333267862 0.005999882151 -0.03933256077 4.357235348e-20 -1.191904678e-10 
1.389595584e-10 1.532416622e-11 -1.786585296e-11 -7.399012931e-11 
-6.346391873e-11
-0.004666575006 -0.01133311073 -0.03799925362 3.09230726e-20 -1.032431091e-10 
1.134837245e-10 -2.140367468e-11 2.352669096e-11 -5.905399583e-11 
-5.372504436e-11

------------without--------------

0.003333267862 0.005999882151 -0.03933256077 2.462301089e+295 2.600734591e+147 
-3.032095915e+147 -3.343731249e+146 3.898326993e+146 -2.200127963e+147 
-1.887126615e+147
-0.004666575006 -0.01133311073 -0.03799925362 inf 4.654777487e+174 
-5.116481775e+174 9.649975086e+173 -1.060714971e+174 -3.75981762e+174 
-3.420536842e+174

Original comment by t.ostase...@gmail.com on 24 Feb 2014 at 5:38

[GoogleCodeExporter]

Resolved,

In bessel function calculations I used absolute value function abs(x), which 
returned int. Changed that now to fabs(x), everything works as expected.

Original comment by t.ostase...@gmail.com on 24 Feb 2014 at 6:19

[GoogleCodeExporter]

The change for this is done now.

Original comment by t.ostase...@gmail.com on 24 Feb 2014 at 6:21

• Changed state: Fixed