Closed GoogleCodeExporter closed 9 years ago
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
> 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
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
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
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
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
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
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
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
The change for this is done now.
Original comment by t.ostase...@gmail.com
on 24 Feb 2014 at 6:21
/cc @to266
Original issue reported on code.google.com by
t.ostase...@gmail.com
on 22 Jan 2014 at 5:59