michael-hartmann / caps Goto Github PK

The integrals for TE and TM polarization are calculated independently at the moment. However, the Legendre polynomials are probably mostly evaluated at the same points. So, caching the values of the associated Legendre polynomials should speed up the calculation drastically.

Update thesis or add new notes that describe what the code is calculating.

http://herbie.uwplse.org

For large values of the imaginary frequency xi, the integrand is basically 0. However, the actual computation may cause problems. So, estimate the value of log(det(Id-M(xi)) for xi>>1 and return 0 for sufficiently high values of xi.

Automatically determine reasonable value of nLeaf

In the MPI program casimir there seems to be some kind of race condition. If there are processes, some values of logdetD(xi,m) for some m are wrong. If only one core is used for the calculation, this problem does not seem to occur.

In section 3.1.1 of the docs it is said that for N cores one should set -n to N+1. The following example contains the option -n 8 which would mean that 7 cores are assumed. This seems a bit odd.

Fix bug. Maybe a problem due to loss of significance?

# compiler: gcc
# compile time: Aug 28 2017 13:55:26
# compiled on: Linux jonas 3.16.0-4-amd64 #1 SMP Debian 3.16.43-2+deb8u3 (2017-08-15) x86_64 GNU/Linux
# LAPACK support: false
# git HEAD: 6fe688a01c5d67de478b7145f8e0eb9bdc864056
# git branch: development
#
# L      = 8e-07
# R      = 0.0001
# LbyR   = 0.008
# T      = 0
# cutoff = 1e-09
# epsrel = 1e-06
# ldim   = 876
# cores  = 8
# alpha  = 0.0158730158730159
# quad   = adaptive Gauss-Kronrod
#
# xi=63, logdetD=-12.0255186746885, t=43.9201
# xi=14683.1056275, logdetD=-6.49554364151339e-101, t=5.1314
# xi=0.270310661838, logdetD=-35.1937568365851, t=44.2067
# xi=2412.82690749, logdetD=-6.68805412840241e-16, t=10.2961
# xi=1.64495844591, logdetD=-34.6906558482894, t=43.6213
# q, a[q], value[q], sign
0 1 314.7357706480661 1
1 -29.40000000000004 317.2108229297982 -1
2 416.2816466552254 318.9527503978982 1
3 -3779.342515860986 320.2476917328405 -1
4 24711.09765451089 321.211764525819 1
5 -123905.1473064856 321.9077043292783 -1
6 495354.2086296655 322.3743744246289 1
7 -1620721.695934563 322.6378096230991 -1
8 4420146.922778942 322.7162899702741 1
9 -10182481.50050869 322.6229902855284 -1
10 20005537.29956692 322.367487405191 1
11 -33755882.5907216 321.9566673299207 -1
12 49153089.52470427 321.3952880883288 1
13 -61955751.84863224 320.6863284213873 -1
14 67699507.73525879 319.831191181788 1
15 -64130282.00927106 318.8297972913867 -1
16 52586264.43578792 317.6805856973787 1
17 -37212346.38851937 316.3804193271365 -1
18 22615661.04308021 314.9243816009094 1
19 -11722423.48194183 313.3054276590975 -1
20 5132415.373043145 311.5138214098096 1
21 -1873142.622510888 309.5362282921238 -1
22 559495.176309977 307.3542079464329 1
23 -133248.6661925225 304.9415647425185 -1
24 24337.12194024635 302.2592853703451 1
25 -3206.758665012491 299.2461392998197 -1
26 445.6078525403483 296.2817189403028 1
27 -17818.24355906366 298.9748709127845 -1
Fatal error: l1=293, l2=27, m=27, p=1, sign=-1, log_I=205.281, q=28 (in _casimir_integrate_I, integration.c:560)

Fatal error: LbyR=0.0008, xi=6102.637461, m=173, ldim=10000 (in slave, casimir_T0.c:412)

The calculations can be sped up drastically. The dominant matrix elements for small separations L/R are near the diagonal.

This has two consequences:

1. long double numerics is no longer needed
2. not all matrix elements have to be calculated

Please consider choosing the latest version of your preferred license.

Get the tests working again.

The function casimir_kernel_M does not handle the conversion from matrix indices i,j to angular momentum numbers l1,l2 dependent on m correctly.

From the documentation it is not really clear what the status of the scripts in the materials directory is. For example, is eps.py supposed to be used by the user and should it be documented then?

I am also wondering about the format of gold.csv. The numerical entries are separated by tabs, so that it is a valid CSV file. However, the documentation says that the entries should be separated by spaces. Reading the CSV file with spaces as separators will not yield the correct result.

I have reread the manuscript. Everything seems fine. The only question is whether the address of @michael-hartmann should be updated in one way or the other. From my side, it would also be fine to leave it as it is.

Using [1] it is possible to calculate ordinary Legendre polynomials in O(1) for arbitrary degree. Using an upwards recurrence relations any Plm can be started from Pl. Surprisingly, the recurrence relations are stable.

Fatal error: ier=2, l=390, m=390, xi_=41943.9, epsrel=4e-09 (in _integrate_K, casimir_scalar.c:82)

The Barnes-Hut algorithm [1], fast multipole methods and hierarchical matrices are somehow related. Try to better understand the connection between those three concepts/ideas. [2] may be useful.

References:
[1] https://www.nature.com/nature/journal/v324/n6096/pdf/324446a0.pdf
[2] Brunner, Junge, Rapp, Bebendorf, Gaul, Comparison of the Fast Multipole Method with Hierarchical Matrices for the Helmholtz-BEM, 2010

As discussed on Firday it might be a good idea to rename this package, because libcasimir is a very general name.

I propose CaPS, Casimir in the Plane-Sphere geometry.

While this name is not very original, it is short, memorable, and catches the essence of the functionality.

Ideas and opinions?

630 /* we could do both calculations together. but it doesn't cost much time -
631 * so why bother?
632 */
633 bessel_lnInuKnu(l-1, chi, &lnIlm, &lnKlm);
634 bessel_lnInuKnu(l, chi, &lnIlp, &lnKlp);

According to the docs and the available flags, it seems that the material properties of the two objects need to be the same. In fact, I added a comment to the docs saying this. But looking at AuAl.c in the examples directory, it seems that different materials can be used. As far as I understand, there is no way to achieve this through flags. Is this correct?

When libblasand liblapack are installed but not libopenblas, the build works fine with the following partial output of CMake:

blas libraries:    /usr/lib/x86_64-linux-gnu/libblas.so
lapack libraries:  /usr/lib/x86_64-linux-gnu/liblapack.so/usr/lib/x86_64-linux-gnu/libblas.so

However, if libopenblas is installed in addition, CMake finds:

blas libraries:    /usr/lib/x86_64-linux-gnu/libopenblas.so
lapack libraries:  /usr/lib/x86_64-linux-gnu/libopenblas.so/usr/lib/x86_64-linux-gnu/libopenblas.so

and the build fails with:

/home/gert/wd/caps/bin/libcaps.so: Warnung: undefinierter Verweis auf »dstemr_«
/home/gert/wd/caps/bin/libcaps.so: Warnung: undefinierter Verweis auf »dpotrf_«
/home/gert/wd/caps/bin/libcaps.so: Warnung: undefinierter Verweis auf »dgemm_«
/home/gert/wd/caps/bin/libcaps.so: Warnung: undefinierter Verweis auf »dgetrf_«
/home/gert/wd/caps/bin/libcaps.so: Warnung: undefinierter Verweis auf »dgeqrf_«
/home/gert/wd/caps/bin/libcaps.so: Warnung: undefinierter Verweis auf »ddot_«
collect2: error: ld returned 1 exit status
CMakeFiles/caps_frontend.dir/build.make:101: recipe for target 'caps' failed
make[2]: *** [caps] Error 1
CMakeFiles/Makefile2:237: recipe for target 'CMakeFiles/caps_frontend.dir/all' failed
make[1]: *** [CMakeFiles/caps_frontend.dir/all] Error 2
Makefile:83: recipe for target 'all' failed
make: *** [all] Error 2

• Write mail to author of HODLR library.
• check and merge fixes.

Try to show that the scattering matrices are positive definite.

First, show that the Mie coefficients al, bl differ in sign. This can be proven using the recurrence relations for PC and probably similar for arbitrary metals.

Then, show that the blocks EE and MM are symmetric and EM = ME and symmetric.

If Id-EE and Id-MM are diagonal dominant, the scattering matrices are positive definite and we can use Choleky decomposition to calculate the determinant.

See Wikipedia article in positive definite matrix.

I tried to follow the compilation procedure described in section 2.1 of the docs and stumbled across a couple of problems (ultimately I succeeded, though).

• cd libcasimir-src/: Probably the directory should be libcasimir-dev aka caps.
• cmake ../src: Since CMakeLists.txt is in libcasimir-dev or whatever is its actual name, this should be cmake .. if called from the subdirectory bin.
• ./casimir_tests: Probably this should be ../caps_tests.

• For submission to JOSS, we need to put the material into a public repo. We could either make this repo public, thereby making the whole project history public or we could create a new repo and transfer the present state there.

• Check numerics (work in progress)
• reimplement the integration for Plasma, xi=0

Use cache.

There may be a problem...

Current status

Perfect mirrors and (partly) Drude mirrors are supported.

Problem

There is no adaptive procedure for integration and no error estimation. There
is no support for real mirrors.

Goal

Add a integration with error estimation, support real mirrors.

Implementation

The Fresnel and Mie coefficients have to be replaced by more general forms. The user should be able to specify an arbitrary dielectric function epsilon(omega).

The integration should be implemented using Gauß-Kronrod. The integral can be split in several subintervals and the values of the associated Legendre polynomials at that nodes can be cached. If the accuracy is not sufficient, the interval with the largest error can be split in two or more intervals to increase accuracy. This process can be repeated until the desired accuracy is reached.

While regenerating output examples, I just noticed that in the high-temperature calculation the start time is recorded in the output but not the stop time. May be the stop time should be added for consistency.

Implement the boundary conditions Y=N and Y=R for the scalar case.

implement a resume feature for casimir

At the moment, the integrals K and I are saved. Each item needs 16 bytes at the moment: 8 bytes double and a 1 byte signed char that sa1610612741ves the sign. With aliasing, this gives 16 bytes. However, the sign is known implicitly, so we don't have to compute it.

Also, maybe we exceed the number of good prime numbers in hash-table.c

Also, we probably don't need to cache every value of I.

Make it possible to select the model for the reflectivity of the mirrors at xi=0. At moment Drude is used.

There seems to be a problem with HODLR:

$ ./caps_logdetD --xi 258.79656952 -L 0.0006481 -R 1 -m 120 --iepsrel 1e-9 --ldim 15430
# ./caps_logdetD, --xi, 258.79656952, -L, 0.0006481, -R, 1, -m, 120, --iepsrel, 1e-9, --ldim, 15430
# L/R    = 0.0006481
# L      = 0.0006481
# R      = 1
# ldim   = 15430
# epsrel = 1.0e-09
# detalg = HODLR
#
# L, R, ξ*(L+R)/c, m, logdet(Id-M), ldim, time
0.0006481, 1, 258.797, 120, -0.004755083928316058, 15430, 53.6748

In contrast with Cholesky logdet=-0.004754880534898247.

Update the HODLR version to the latest available on github. The latest version on github includes a version for symmetric, positive definite matrices.

Also, update libeigen to the latest version.

Create data for Drude metals at T=0. (Drude-Model, not optical data)

Rereading the manual, I found the following open issues which should be addressed in one way or the other before submission:

• The manual contains different version numbers (0.4.2 on the title page and 0.4.3 in the listings). Before fixing this, we should decide whether in view of the submission, we want to go to 1.0. (resolved by #123 and #126, sticking to 0.4.3)
• Is the reference to Emig et al. in the context of plane-cylinder really the appropriate one? (fixed by #124)
• Is it ok to leave local computer information in the manual? (addressed in #127)
• The syntax highlighting in some cases seems a bit odd. For example, on page 5, the argument of -n is highlighted but not the other numerical parameters. Overall, the highlighting seems pretty inconsistent and in most cases rather useless. Should we get rid of the highlighting or rather keep it? (fixed by #125)

I can do the changes once it is decided what to change or keep.

In order to better understand low-rank approximations, look at the rank of different sub-blocks, and have a look at the tree that HODLR creates (header/HODLR_Node.hpp).

In the numerics the signs can be drastically simplified for perfect reflectors. It is though an open question if this is also true for Drude metals and materials with a generic dielectric function.

One can show that
l_I_{l+1/2}(chi) - chi_I_{l-1/2}(chi) < 0
for all chi and l. Therefore, for perfect reflectors the sign of
al is give by (-1)^l and the sign of b_l is given by (-1)^(l+1}.

One can also show that:
sla < 0
slb < 0
slc > 0
sld < 0

Therefore th denominators of the Mie coefficients,
slc-sld
n^2*slc-sld
are always positive. However, it is not clear what the sign of
the nominators of a_l and b_l are for arbitrary epsilon.

https://refractiveindex.info/

The matrices are banded and we don't exploit this property at the moment. Using dgbtrf would speed up calculation.

See http://physics.oregonstate.edu/~landaur/nacphy/lapack/routines/dgbtrf.html