From 20fac6036fa8e0a2d7390bfdae685d9e06852e4e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 19 Aug 2019 15:10:39 +0300 Subject: [PATCH] Complex wave number for ewald sum (failing) + unit test. Former-commit-id: a6c5e566028602ecd1047bc0b44677a4eeae3046 --- qpms/ewald.c | 4 +- qpms/translations.c | 3 +- qpms/translations.h | 33 +---- tests/CMakeLists.txt | 14 ++ tests/test_scalar_ewald32.c | 283 ++++++++++++++++++++++++++++++++++++ 5 files changed, 301 insertions(+), 36 deletions(-) create mode 100644 tests/test_scalar_ewald32.c diff --git a/qpms/ewald.c b/qpms/ewald.c index 2bde7b9..26bda5d 100644 --- a/qpms/ewald.c +++ b/qpms/ewald.c @@ -6,6 +6,7 @@ #include #include #include "tiny_inlines.h" +#include "qpms_error.h" #include #include #include @@ -173,7 +174,6 @@ void qpms_ewald3_constants_free(qpms_ewald3_constants_t *c) { } - int ewald3_sigma0(complex double *result, double *err, const qpms_ewald3_constants_t *c, const double eta, const complex double k) @@ -467,7 +467,7 @@ int ewald3_sigma_long ( return ewald3_1_z_sigma_long(target, err, c, eta, k, unitcell_volume, latdim, pgen_K, pgen_generates_shifted_points, beta, particle_shift); // TODO 3D case and general 2D case - else abort(); // NOT IMPLEMENTED + else QPMS_NOT_IMPLEMENTED("3D or general 2D (outside XY plane) Ewald sum."); } struct sigma2_integrand_params { diff --git a/qpms/translations.c b/qpms/translations.c index 062aff4..2832cdc 100644 --- a/qpms/translations.c +++ b/qpms/translations.c @@ -857,7 +857,6 @@ int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_tr #ifdef LATTICESUMS32 - // N.B. alternative point generation strategy toggled by macro GEN_RSHIFTEDPOINTS // and GEN_KSHIFTEDPOINTS. // The results should be the same. The performance can slightly differ (especially @@ -867,7 +866,7 @@ int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c, complex double * const Bdest, double * const Berr, const ptrdiff_t deststride, const ptrdiff_t srcstride, /* qpms_bessel_t J*/ // assume QPMS_HANKEL_PLUS - const double eta, const double k, + const double eta, const complex double k, const cart2_t b1, const cart2_t b2, const cart2_t beta, const cart2_t particle_shift, diff --git a/qpms/translations.h b/qpms/translations.h index 55c4bd3..ef8cb48 100644 --- a/qpms/translations.h +++ b/qpms/translations.h @@ -154,43 +154,12 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p( // TODO more high-level functions with more advanced lattice generators etc. (after // the prerequisities from lattices2d.h are implememted) -#if 0 // NI -int qpms_trans_calculator_e32_long_points_and_shift(const qpms_trans_calculator *c, - complex double *Adest_long, double *Aerr_long, - complex double *Bdest_long, double *Berr_long, - double eta, double k, double unitcell_area, - size_t npoints, const cart2_t *Kpoints, - cart2_t beta, - cart2_t particle_shift - ); - -int qpms_trans_calculator_e32_short_points_and_shift(const qpms_trans_calculator *c, - complex double *Adest_short, double *Aerr_short, - complex double *Bdest_short, double *Berr_short, - double eta, double k, - size_t npoints, const cart2_t *Rpoints, - cart2_t beta, - cart2_t particle_shift - ); -#endif - -int qpms_trans_calculator_get_AB_arrays_e32_both_points_and_shift(const qpms_trans_calculator *c, - complex double *Adest, double *Aerr, - complex double *Bdest, double *Berr, - const ptrdiff_t deststride, const ptrdiff_t srcstride, - const double eta, const double k, - const double unitcell_area, - const size_t nRpoints, const cart2_t *Rpoints, - const size_t nKpoints, const cart2_t *Kpoints, - const cart2_t beta, - const cart2_t particle_shift - ); int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c, complex double *Adest, double *Aerr, complex double *Bdest, double *Berr, const ptrdiff_t deststride, const ptrdiff_t srcstride, - const double eta, const double k, + const double eta, const complex double k, cart2_t b1, cart2_t b2, const cart2_t beta, const cart2_t particle_shift, diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 61ae546..e42300b 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -11,6 +11,20 @@ add_executable( test_single_translations_vs_calc single_translations_vs_calc.c ) target_link_libraries( test_single_translations_vs_calc qpms gsl lapacke amos m ) target_include_directories( test_single_translations_vs_calc PRIVATE .. ) +add_executable( test_scalar_ewald32 test_scalar_ewald32.c ) +target_link_libraries( test_scalar_ewald32 qpms gsl lapacke amos m ) +target_include_directories( test_scalar_ewald32 PRIVATE .. ) + add_custom_target( mytests DEPENDS test_single_translations_vs_calc test_vswf_translations test_vswf_translations_array ) add_test( NAME single_vs_array_translation_coeffs COMMAND test_single_translations_vs_calc ) + +add_test( NAME scalar_ewald32_realk1 COMMAND test_scalar_ewald32 +# lMax b1.x b1.y b2.x b2.y wavenum.real wavenum.imag k.x k.y particle_shift.x particle_shift.y csphase rtol atol maxR maxK eta1 ... + 3 1 0 0 1 2.3 0 2.7 1 0.5 0.1325 -1 1e-8 1e-10 20 160 1.5 1.6 2.5 2.6 + ) + +add_test( NAME scalar_ewald32_cplxk1 COMMAND test_scalar_ewald32 +# lMax b1.x b1.y b2.x b2.y wavenum.real wavenum.imag k.x k.y particle_shift.x particle_shift.y csphase rtol atol maxR maxK eta1 ... + 3 1 0 0 1 2.3 0.1 2.7 1 0.5 0.1325 -1 1e-8 1e-10 20 160 1.5 1.6 2.5 2.6 + ) diff --git a/tests/test_scalar_ewald32.c b/tests/test_scalar_ewald32.c new file mode 100644 index 0000000..16a9291 --- /dev/null +++ b/tests/test_scalar_ewald32.c @@ -0,0 +1,283 @@ +// Perform Ewald summation (2D xy-lattice in 3D space ) of SSWFs with different Ewald parameters and check whether the difference is inside the tolerance range. +// run as +// test_scalar_ewald32 lMax a1.x a1.y a2.x a2.y wavenum.real wavenum.imag k.x k.y particle_shift.x particle_shift.y csphase rtol atol maxR maxK eta1 [eta2 [eta 3 ...]] +// c99 -o ewaldshift3g_vargeom -ggdb -Wall -I ../ ewaldshift3g_vargeom.c ../qpms/ewald.c ../qpms/ewaldsf.c ../qpms/lattices2d.c ../qpms/latticegens.c -lgsl -lm -lblas + + + +// implementation of the [LT(4.16)] test +#include +#define M_SQRTPI 1.7724538509055160272981674833411452 +#define M_SQRT3 1.7320508075688772935274463415058724 +#include +#include +#include +#include +#include +#include +#include +typedef struct ewaldtest2d_params { + qpms_l_t lMax; + point2d b1, b2; + point2d beta; + point2d particle_shift; + complex double k; + //double a; + double eta; + double maxR; + double maxK; + double csphase; +} ewaldtest2d_params; + +typedef struct ewaldtest2d_results { + ewaldtest2d_params p; + complex double *sigmas_short, + *sigmas_long, + sigma0, + *sigmas_total; + double *err_sigmas_short, + *err_sigmas_long, + err_sigma0, + *err_sigmas_total; + complex double *regsigmas_416; +} ewaldtest2d_results; + + +void ewaldtest2d_results_free(ewaldtest2d_results *r) { + free(r->sigmas_short); + free(r->sigmas_long); + free(r->sigmas_total); + free(r->err_sigmas_long); + free(r->err_sigmas_total); + free(r->err_sigmas_short); + free(r->regsigmas_416); + free(r); +} + +static inline double san(double x) { + return fabs(x) < 1e-13 ? 0 : x; +} + +int isclose_cmplx(complex double a, complex double b, double rtol, double atol) { + return cabs(a-b) <= atol + rtol * .5 * (cabs(a) + cabs(b)); +} + +ewaldtest2d_results *ewaldtest2d(const ewaldtest2d_params p); + +int main(int argc, char **argv) { + bool verbose = !!getenv("QPMS_VERBOSE_TESTS"); + gsl_set_error_handler(IgnoreUnderflowsGSLErrorHandler); + QPMS_ENSURE(argc >= 18, "At least 16 arguments expected, I see only %d.", argc-1); + int netas = argc - 17; + ewaldtest2d_params plist[netas]; + double atol, rtol; + plist[0].lMax = atoi(argv[1]); + plist[0].b1.x = strtod(argv[2], NULL); + plist[0].b1.y = strtod(argv[3], NULL); + plist[0].b2.x = strtod(argv[4], NULL); + plist[0].b2.y = strtod(argv[5], NULL); + plist[0].k = strtod(argv[6], NULL) + I*strtod(argv[7], NULL); + plist[0].beta.x = strtod(argv[8], NULL); + plist[0].beta.y = strtod(argv[9], NULL); + plist[0].particle_shift.x = strtod(argv[10], NULL); + plist[0].particle_shift.y = strtod(argv[11], NULL); + plist[0].csphase = strtod(argv[12], NULL); + atol = strtod(argv[13], NULL); + rtol = strtod(argv[14], NULL); + plist[0].maxR = strtod(argv[15], NULL); + plist[0].maxK = strtod(argv[16], NULL); + plist[0].eta = strtod(argv[17], NULL); + for(int i = 1; i < netas; ++i) { + plist[i] = plist[0]; + plist[i].eta = strtod(argv[17+i], NULL); + } + + ewaldtest2d_results *r[netas]; + + int fails = 0; + + for (size_t i = 0; i < netas; ++i) { + ewaldtest2d_params p = plist[i]; + r[i] = ewaldtest2d(p); + // TODO print per-test header here + printf("===============================\n"); + printf("b1 = (%g, %g), b2 = (%g, %g)," /* "K1 = (%g, %g), K2 = (%g, %g),"*/ " Kmax = %g, Rmax = %g, lMax = %d, eta = %g, k = %g%+gj, beta = (%g,%g), ps = (%g,%g), csphase = %g\n", + p.b1.x, p.b1.y, p.b2.x, p.b2.y,/*TODO K1, K2*/ p.maxK, p.maxR, p.lMax, p.eta, creal(p.k), cimag(p.k), p.beta.x, p.beta.y, p.particle_shift.x, p.particle_shift.y, p.csphase); + printf("sigma0: %.16g%+.16gj\n", creal(r[i]->sigma0), cimag(r[i]->sigma0)); + for (qpms_l_t n = 0; n <= p.lMax; ++n) { + for (qpms_m_t m = -n; m <= n; ++m){ + if ((m+n)%2) continue; + qpms_y_t y = qpms_mn2y_sc(m,n); + qpms_y_t y_conj = qpms_mn2y_sc(-m,n); + // y n m sigma_total (err), regsigmas_416 regsigmas_415_recon + if (verbose) printf("%zd %d %d: T:%.16g%+.16gj(%.3g) L:%.16g%+.16gj(%.3g) S:%.16g%+.16gj(%.3g) \n" + //"| predict %.16g%+.16gj \n| actual %.16g%+.16gj\n" + , + y, n, m, creal(san(r[i]->sigmas_total[y])), san(cimag(r[i]->sigmas_total[y])), + r[i]->err_sigmas_total[y], + san(creal(r[i]->sigmas_long[y])), san(cimag(r[i]->sigmas_long[y])), + r[i]->err_sigmas_long[y], + san(creal(r[i]->sigmas_short[y])), san(cimag(r[i]->sigmas_short[y])), + r[i]->err_sigmas_short[y] + // TODO and count big differences as failures. + //san(creal(r[i]->regsigmas_416[y])), san(cimag(r[i]->regsigmas_416[y])), + //san(creal(r[i]->sigmas_total[y]) + creal(r[i]->sigmas_total[y_conj])), + //san(cimag(r[i]->sigmas_total[y]) - cimag(r[i]->sigmas_total[y_conj])) + ); + } + } + } + for (size_t i = 0; i < netas; ++i) { + for (size_t j = i+1; j < netas; ++j){ + for (qpms_l_t n = 0; n <= plist[i].lMax; ++n) { + for (qpms_m_t m = -n; m <= n; ++m){ + if ((m+n)%2) continue; + qpms_y_t y = qpms_mn2y_sc(m,n); + qpms_y_t y_conj = qpms_mn2y_sc(-m,n); + if (!isclose_cmplx(r[i]->sigmas_total[y], r[j]->sigmas_total[y], rtol, atol)) { + ++fails; + printf("with eta = %.16g:\n", plist[i].eta); + printf("%zd %d %d: T:%.16g%+.16gj(%.3g) L:%.16g%+.16gj(%.3g) S:%.16g%+.16gj(%.3g) \n" + //"| predict %.16g%+.16gj \n| actual %.16g%+.16gj\n" + , + y, n, m, creal(san(r[i]->sigmas_total[y])), san(cimag(r[i]->sigmas_total[y])), + r[i]->err_sigmas_total[y], + san(creal(r[i]->sigmas_long[y])), san(cimag(r[i]->sigmas_long[y])), + r[i]->err_sigmas_long[y], + san(creal(r[i]->sigmas_short[y])), san(cimag(r[i]->sigmas_short[y])), + r[i]->err_sigmas_short[y] + // TODO and count big differences as failures. + //san(creal(r[i]->regsigmas_416[y])), san(cimag(r[i]->regsigmas_416[y])), + //san(creal(r[i]->sigmas_total[y]) + creal(r[i]->sigmas_total[y_conj])), + //san(cimag(r[i]->sigmas_total[y]) - cimag(r[i]->sigmas_total[y_conj])) + ); + printf("with eta = %.16g:\n", plist[j].eta); + printf("%zd %d %d: T:%.16g%+.16gj(%.3g) L:%.16g%+.16gj(%.3g) S:%.16g%+.16gj(%.3g) \n" + //"| predict %.16g%+.16gj \n| actual %.16g%+.16gj\n" + , + y, n, m, creal(san(r[j]->sigmas_total[y])), san(cimag(r[j]->sigmas_total[y])), + r[j]->err_sigmas_total[y], + san(creal(r[j]->sigmas_long[y])), san(cimag(r[j]->sigmas_long[y])), + r[j]->err_sigmas_long[y], + san(creal(r[j]->sigmas_short[y])), san(cimag(r[j]->sigmas_short[y])), + r[j]->err_sigmas_short[y] + // TODO and count big differences as failures. + //san(creal(r[j]->regsigmas_416[y])), san(cimag(r[j]->regsigmas_416[y])), + //san(creal(r[j]->sigmas_total[y]) + creal(r[j]->sigmas_total[y_conj])), + //san(cimag(r[j]->sigmas_total[y]) - cimag(r[j]->sigmas_total[y_conj])) + ); + } + } + } + } + ewaldtest2d_results_free(r[i]); + } + + return fails; +} + + +int ewaldtest_counter = 0; + + +ewaldtest2d_results *ewaldtest2d(const ewaldtest2d_params p) { + cart3_t beta3 = cart22cart3xy(p.beta); + cart3_t particle_shift3 = cart22cart3xy(p.particle_shift); + + cart2_t b1 = p.b1, b2 = p.b2, rb1, rb2; + if (QPMS_SUCCESS != l2d_reciprocalBasis2pi(b1, b2, &rb1, &rb2)) + abort(); + + const double A = l2d_unitcell_area(b1, b2); // sqrt(3) * a * a / 2.; // unit cell size + const double K_len = cart2norm(rb1)+cart2norm(rb2); //4*M_PI/a/sqrt(3); // reciprocal vector length + + ewaldtest2d_results *results = malloc(sizeof(ewaldtest2d_results)); + results->p = p; + + // skip zeroth point if it coincides with origin + bool include_origin = !(fabs(p.particle_shift.x) == 0 + && fabs(p.particle_shift.y) == 0); + + PGen Rlgen = PGen_xyWeb_new(b1, b2, BASIS_RTOL, CART2_ZERO, 0, include_origin, p.maxR, false); + //PGen Rlgen_plus_shift = PGen_xyWeb_new(b1, b2, BASIS_RTOL, cart2_scale(-1 /* CHECKSIGN */, particle_shift2), 0, include_origin, p.maxR + a, false); + PGen Klgen = PGen_xyWeb_new(rb1, rb2, BASIS_RTOL, CART2_ZERO, 0, true, p.maxK + K_len, false); + //PGen Klgen_plus_beta = PGen_xyWeb_new(rb1, rb2, BASIS_RTOL, beta2, 0, true, p.maxK + K_len, false); + + qpms_y_t nelem_sc = qpms_lMax2nelem_sc(p.lMax); + + results->sigmas_short = malloc(sizeof(complex double)*nelem_sc); + results->sigmas_long = malloc(sizeof(complex double)*nelem_sc); + results->sigmas_total = malloc(sizeof(complex double)*nelem_sc); + results->err_sigmas_short = malloc(sizeof(double)*nelem_sc); + results->err_sigmas_long = malloc(sizeof(double)*nelem_sc); + results->err_sigmas_total = malloc(sizeof(double)*nelem_sc); + + qpms_ewald3_constants_t *c = qpms_ewald3_constants_init(p.lMax, p.csphase); + + if (0!=ewald3_sigma_long(results->sigmas_long, + results->err_sigmas_long, c, p.eta, p.k, A, + LAT_2D_IN_3D_XYONLY, &Klgen, false, beta3, particle_shift3)) + abort(); + if (0!=ewald3_sigma_short( + results->sigmas_short, results->err_sigmas_short, c, + p.eta, p.k, LAT_2D_IN_3D_XYONLY, &Rlgen, false, beta3, particle_shift3)) + abort(); + if (0!=ewald3_sigma0(&(results->sigma0), &(results->err_sigma0), c, p.eta, p.k)) + abort(); + for(qpms_y_t y = 0; y < nelem_sc; ++y) { + results->sigmas_total[y] = results->sigmas_short[y] + results->sigmas_long[y]; + results->err_sigmas_total[y] = results->err_sigmas_short[y] + results->err_sigmas_long[y]; + } + results->sigmas_total[0] += results->sigma0; + results->err_sigmas_total[0] += results->err_sigma0; + + // Now calculate the reference values [LT(4.16)] + results->regsigmas_416 = calloc(nelem_sc, sizeof(complex double)); + results->regsigmas_416[0] = -2 * c->legendre0[gsl_sf_legendre_array_index(0,0)]; + +#if 0 // not yet implemented for the new API + { + double legendres[gsl_sf_legendre_array_n(p.lMax)]; + points2d_rordered_t sel = + points2d_rordered_annulus(Kpoints_plus_beta, 0, true, p.k, false); + if (0 != sel.nrs) + { + point2d *beta_pq_lessthan_k = sel.base + sel.r_offsets[0]; + size_t beta_pq_lessthan_k_count = sel.r_offsets[sel.nrs] - sel.r_offsets[0]; + for(size_t i = 0; i < beta_pq_lessthan_k_count; ++i) { + point2d beta_pq = beta_pq_lessthan_k[i]; + double rbeta_pq = cart2norm(beta_pq); + double arg_pq = atan2(beta_pq.y, beta_pq.x); + double denom = sqrt(p.k*p.k - rbeta_pq*rbeta_pq); + if( gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE, + p.lMax, denom/p.k, p.csphase, legendres) != 0) + abort(); + for (qpms_y_t y = 0; y < nelem_sc; ++y) { + qpms_l_t n; qpms_m_t m; + qpms_y2mn_sc_p(y, &m, &n); + if ((m+n)%2 != 0) + continue; + complex double eimf = cexp(I*m*arg_pq); + results->regsigmas_416[y] += + 4*M_PI*ipow(n)/p.k/A + * eimf * legendres[gsl_sf_legendre_array_index(n,abs(m))] * min1pow_m_neg(m) + / denom; + } + } + } + } +#else + for(qpms_y_t y = 0; y < nelem_sc; ++y) { + qpms_l_t n; qpms_m_t m; + qpms_y2mn_sc_p(y, &m, &n); + if ((m+n)%2 != 0) + continue; + results->regsigmas_416[y] = NAN; + } +#endif + + + qpms_ewald3_constants_free(c); + ++ewaldtest_counter; + return results; +}