tests/test_scalar_ewald32.c

   1 // 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.
   2 // run as
   3 // 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 ...]]
   4 // 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
   5
   6
   7
   8 // implementation of the [LT(4.16)] test
   9 #include <math.h>
  10 #define M_SQRTPI 1.7724538509055160272981674833411452
  11 #define M_SQRT3 1.7320508075688772935274463415058724
  12 #include <qpms/ewald.h>
  13 #include <qpms/tiny_inlines.h>
  14 #include <qpms/indexing.h>
  15 #include <stdlib.h>
  16 #include <stdio.h>
  17 #include <string.h>
  18 #include <float.h>
  19 #include <gsl/gsl_sf_legendre.h>
  20 #include <fenv.h>
  21
  22 typedef struct ewaldtest2d_params {
  23   qpms_l_t lMax;
  24   point2d b1, b2;
  25   point2d beta;
  26   point2d particle_shift;
  27   complex double k;
  28   //double a;
  29   double eta;
  30   double maxR;
  31   double maxK;
  32   double csphase;
  33 } ewaldtest2d_params;
  34
  35 typedef struct ewaldtest2d_results {
  36   ewaldtest2d_params p;
  37   complex double *sigmas_short,
  38           *sigmas_long,
  39           sigma0,
  40           *sigmas_total;
  41   double *err_sigmas_short,
  42     *err_sigmas_long,
  43     err_sigma0,
  44     *err_sigmas_total;
  45   complex double *regsigmas_416;
  46 } ewaldtest2d_results;
  47
  48
  49 void ewaldtest2d_results_free(ewaldtest2d_results *r) {
  50   free(r->sigmas_short);
  51   free(r->sigmas_long);
  52   free(r->sigmas_total);
  53   free(r->err_sigmas_long);
  54   free(r->err_sigmas_total);
  55   free(r->err_sigmas_short);
  56   free(r->regsigmas_416);
  57   free(r);
  58 }
  59
  60 static inline double san(double x) {
  61   return fabs(x) < 1e-13 ? 0 : x;
  62 }
  63
  64 int isclose_cmplx(complex double a, complex double b, double rtol, double atol) {
  65   return cabs(a-b) <= atol + rtol * .5 * (cabs(a) + cabs(b));
  66 }
  67
  68 ewaldtest2d_results *ewaldtest2d(const ewaldtest2d_params p);
  69
  70 int main(int argc, char **argv) {
  71   feenableexcept(FE_INVALID | FE_OVERFLOW);
  72   bool verbose = !!getenv("QPMS_VERBOSE_TESTS");
  73   gsl_set_error_handler(IgnoreUnderflowsGSLErrorHandler);
  74   QPMS_ENSURE(argc >= 18, "At least 16 arguments expected, I see only %d.", argc-1);
  75   int netas = argc - 17;
  76   ewaldtest2d_params plist[netas];
  77   double atol, rtol;
  78   plist[0].lMax = atoi(argv[1]);
  79   plist[0].b1.x = strtod(argv[2], NULL);
  80   plist[0].b1.y = strtod(argv[3], NULL);
  81   plist[0].b2.x = strtod(argv[4], NULL);
  82   plist[0].b2.y = strtod(argv[5], NULL);
  83   plist[0].k = strtod(argv[6], NULL) + I*strtod(argv[7], NULL);
  84   plist[0].beta.x = strtod(argv[8], NULL);
  85   plist[0].beta.y = strtod(argv[9], NULL);
  86   plist[0].particle_shift.x = strtod(argv[10], NULL);
  87   plist[0].particle_shift.y = strtod(argv[11], NULL);
  88   plist[0].csphase = strtod(argv[12], NULL);
  89   atol = strtod(argv[13], NULL);
  90   rtol = strtod(argv[14], NULL);
  91   plist[0].maxR = strtod(argv[15], NULL);
  92   plist[0].maxK = strtod(argv[16], NULL);
  93   plist[0].eta  = strtod(argv[17], NULL);
  94   for(int i = 1; i < netas; ++i) {
  95     plist[i] = plist[0];
  96     plist[i].eta = strtod(argv[17+i], NULL);
  97   }
  98
  99   ewaldtest2d_results *r[netas];
 100
 101   int fails = 0;
 102
 103   for (size_t i = 0; i < netas; ++i) {
 104     ewaldtest2d_params p = plist[i];
 105     r[i] = ewaldtest2d(p);
 106     // TODO print per-test header here
 107     printf("===============================\n");
 108     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",
 109         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);
 110     printf("sigma0: %.16g%+.16gj\n", creal(r[i]->sigma0), cimag(r[i]->sigma0));
 111     for (qpms_l_t n = 0; n <= p.lMax; ++n) {
 112       for (qpms_m_t m = -n; m <= n; ++m){
 113         if ((m+n)%2) continue;
 114         qpms_y_t y = qpms_mn2y_sc(m,n);
 115         qpms_y_t y_conj = qpms_mn2y_sc(-m,n);
 116         // y n m sigma_total (err), regsigmas_416 regsigmas_415_recon
 117         if (verbose) printf("%zd %d %d: T:%.16g%+.16gj(%.3g) L:%.16g%+.16gj(%.3g) S:%.16g%+.16gj(%.3g) \n"
 118             //"| predict %.16g%+.16gj \n| actual  %.16g%+.16gj\n"
 119             ,
 120             y, n, m, creal(san(r[i]->sigmas_total[y])), san(cimag(r[i]->sigmas_total[y])),
 121             r[i]->err_sigmas_total[y],
 122             san(creal(r[i]->sigmas_long[y])), san(cimag(r[i]->sigmas_long[y])),
 123             r[i]->err_sigmas_long[y],
 124             san(creal(r[i]->sigmas_short[y])), san(cimag(r[i]->sigmas_short[y])),
 125             r[i]->err_sigmas_short[y]
 126             // TODO and count big differences as failures.
 127             //san(creal(r[i]->regsigmas_416[y])), san(cimag(r[i]->regsigmas_416[y])),
 128             //san(creal(r[i]->sigmas_total[y]) + creal(r[i]->sigmas_total[y_conj])),
 129             //san(cimag(r[i]->sigmas_total[y]) - cimag(r[i]->sigmas_total[y_conj]))
 130         );
 131       }
 132     }
 133   }
 134
 135   bool toprint[netas];
 136
 137   for (qpms_l_t n = 0; n <= plist[0].lMax; ++n) {
 138     for (qpms_m_t m = -n; m <= n; ++m){
 139       memset(toprint, 0, netas*sizeof(bool));
 140       if ((m+n)%2) continue;
 141       qpms_y_t y = qpms_mn2y_sc(m,n);
 142       qpms_y_t y_conj = qpms_mn2y_sc(-m,n);
 143       for (size_t i = 0; i < netas; ++i) {
 144         for (size_t j = i+1; j < netas; ++j){
 145           if (!isclose_cmplx(r[i]->sigmas_total[y], r[j]->sigmas_total[y], rtol, atol))
 146             toprint[i] = toprint[j] = true;
 147         }
 148         if (toprint[i]) {
 149           ++fails;
 150           printf("with eta = %.16g:\n", plist[i].eta);
 151           printf("%zd %d %d: T:%.16g%+.16gj(%.3g) L:%.16g%+.16gj(%.3g) S:%.16g%+.16gj(%.3g) \n"
 152               //"| predict %.16g%+.16gj \n| actual  %.16g%+.16gj\n"
 153               ,
 154               y, n, m, creal(san(r[i]->sigmas_total[y])), san(cimag(r[i]->sigmas_total[y])),
 155               r[i]->err_sigmas_total[y],
 156               san(creal(r[i]->sigmas_long[y])), san(cimag(r[i]->sigmas_long[y])),
 157               r[i]->err_sigmas_long[y],
 158               san(creal(r[i]->sigmas_short[y])), san(cimag(r[i]->sigmas_short[y])),
 159               r[i]->err_sigmas_short[y]
 160               // TODO and count big differences as failures.
 161               //san(creal(r[i]->regsigmas_416[y])), san(cimag(r[i]->regsigmas_416[y])),
 162               //san(creal(r[i]->sigmas_total[y]) + creal(r[i]->sigmas_total[y_conj])),
 163               //san(cimag(r[i]->sigmas_total[y]) - cimag(r[i]->sigmas_total[y_conj]))
 164               );
 165           //if(!y) printf("0:%.16g%+.16g\n", san(creal(r[i]->sigma0)), san(cimag(r[i]->sigma0)));
 166         }
 167       }
 168     }
 169   }
 170
 171   for (size_t i = 0; i < netas; ++i)
 172     ewaldtest2d_results_free(r[i]);
 173
 174   return fails;
 175 }
 176
 177
 178 int ewaldtest_counter = 0;
 179
 180
 181 ewaldtest2d_results *ewaldtest2d(const ewaldtest2d_params p) {
 182   cart3_t beta3 = cart22cart3xy(p.beta);
 183   cart3_t particle_shift3 = cart22cart3xy(p.particle_shift);
 184
 185   cart2_t b1 = p.b1, b2 = p.b2, rb1, rb2;
 186   if (QPMS_SUCCESS != l2d_reciprocalBasis2pi(b1, b2, &rb1, &rb2))
 187     abort();
 188
 189   const double A = l2d_unitcell_area(b1, b2); // sqrt(3) * a * a / 2.; // unit cell size
 190   const double K_len = cart2norm(rb1)+cart2norm(rb2); //4*M_PI/a/sqrt(3); // reciprocal vector length
 191
 192   ewaldtest2d_results *results = malloc(sizeof(ewaldtest2d_results));
 193   results->p = p;
 194
 195   // skip zeroth point if it coincides with origin
 196   bool include_origin = !(fabs(p.particle_shift.x) == 0
 197       && fabs(p.particle_shift.y) == 0);
 198
 199   PGen Rlgen = PGen_xyWeb_new(b1, b2, BASIS_RTOL, CART2_ZERO, 0, include_origin, p.maxR, false);
 200   //PGen Rlgen_plus_shift = PGen_xyWeb_new(b1, b2, BASIS_RTOL, cart2_scale(-1 /* CHECKSIGN */, particle_shift2), 0, include_origin, p.maxR + a, false);
 201   PGen Klgen = PGen_xyWeb_new(rb1, rb2, BASIS_RTOL, CART2_ZERO, 0, true, p.maxK + K_len, false);
 202   //PGen Klgen_plus_beta = PGen_xyWeb_new(rb1, rb2, BASIS_RTOL, beta2, 0, true, p.maxK + K_len, false);
 203
 204   qpms_y_t nelem_sc = qpms_lMax2nelem_sc(p.lMax);
 205
 206   results->sigmas_short = malloc(sizeof(complex double)*nelem_sc);
 207   results->sigmas_long = malloc(sizeof(complex double)*nelem_sc);
 208   results->sigmas_total = malloc(sizeof(complex double)*nelem_sc);
 209   results->err_sigmas_short = malloc(sizeof(double)*nelem_sc);
 210   results->err_sigmas_long = malloc(sizeof(double)*nelem_sc);
 211   results->err_sigmas_total = malloc(sizeof(double)*nelem_sc);
 212
 213   qpms_ewald3_constants_t *c = qpms_ewald3_constants_init(p.lMax, p.csphase);
 214
 215   if (0!=ewald3_sigma_long(results->sigmas_long,
 216         results->err_sigmas_long, c, p.eta, p.k, A,
 217         LAT_2D_IN_3D_XYONLY, &Klgen, false, beta3, particle_shift3))
 218       abort();
 219   if (0!=ewald3_sigma_short(
 220     results->sigmas_short, results->err_sigmas_short, c,
 221     p.eta, p.k, LAT_2D_IN_3D_XYONLY, &Rlgen, false,  beta3, particle_shift3))
 222     abort();
 223   if (0!=ewald3_sigma0(&(results->sigma0), &(results->err_sigma0), c, p.eta, p.k))
 224     abort();
 225   for(qpms_y_t y = 0; y < nelem_sc; ++y) {
 226     results->sigmas_total[y] = results->sigmas_short[y] + results->sigmas_long[y];
 227     results->err_sigmas_total[y] = results->err_sigmas_short[y] + results->err_sigmas_long[y];
 228   }
 229   if(!include_origin) { // "Renormalised" contribution of origin point
 230     results->sigmas_total[0] += results->sigma0;
 231     results->err_sigmas_total[0] += results->err_sigma0;
 232   }
 233
 234   // Now calculate the reference values [LT(4.16)]
 235   results->regsigmas_416 = calloc(nelem_sc, sizeof(complex double));
 236   results->regsigmas_416[0] = -2 * c->legendre0[gsl_sf_legendre_array_index(0,0)];
 237
 238 #if 0 // not yet implemented for the new API
 239   {
 240     double legendres[gsl_sf_legendre_array_n(p.lMax)];
 241     points2d_rordered_t sel =
 242       points2d_rordered_annulus(Kpoints_plus_beta, 0, true, p.k, false);
 243     if (0 != sel.nrs)
 244     {
 245       point2d *beta_pq_lessthan_k = sel.base + sel.r_offsets[0];
 246       size_t beta_pq_lessthan_k_count = sel.r_offsets[sel.nrs] - sel.r_offsets[0];
 247       for(size_t i = 0; i < beta_pq_lessthan_k_count; ++i) {
 248         point2d beta_pq = beta_pq_lessthan_k[i];
 249         double rbeta_pq = cart2norm(beta_pq);
 250         double arg_pq = atan2(beta_pq.y, beta_pq.x);
 251         double denom = sqrt(p.k*p.k - rbeta_pq*rbeta_pq);
 252         if( gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,
 253               p.lMax, denom/p.k, p.csphase, legendres) != 0)
 254           abort();
 255         for (qpms_y_t y = 0; y < nelem_sc; ++y) {
 256           qpms_l_t n; qpms_m_t m;
 257           qpms_y2mn_sc_p(y, &m, &n);
 258           if ((m+n)%2 != 0)
 259             continue;
 260           complex double eimf = cexp(I*m*arg_pq);
 261           results->regsigmas_416[y] +=
 262             4*M_PI*ipow(n)/p.k/A
 263             * eimf * legendres[gsl_sf_legendre_array_index(n,abs(m))] * min1pow_m_neg(m)
 264             / denom;
 265         }
 266       }
 267     }
 268   }
 269 #else
 270   for(qpms_y_t y = 0; y < nelem_sc; ++y) {
 271           qpms_l_t n; qpms_m_t m;
 272           qpms_y2mn_sc_p(y, &m, &n);
 273           if ((m+n)%2 != 0)
 274             continue;
 275           results->regsigmas_416[y] = NAN;
 276   }
 277 #endif
 278
 279
 280   qpms_ewald3_constants_free(c);
 281   ++ewaldtest_counter;
 282   return results;
 283 }