From d5af94903245873fcf2797ad80648e9cd92f735e Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Sun, 2 Aug 2020 11:05:47 +0200
Subject: [PATCH] fix new symmetrize

---
 TurTLE/test/test_fftw.py          | 36 ++++++++++++++++-
 cpp/field.cpp                     | 17 ++++----
 cpp/field.hpp                     |  2 +-
 cpp/full_code/symmetrize_test.cpp | 67 +++++++++++--------------------
 4 files changed, 68 insertions(+), 54 deletions(-)

diff --git a/TurTLE/test/test_fftw.py b/TurTLE/test/test_fftw.py
index e03d93f9..1b1d4e63 100644
--- a/TurTLE/test/test_fftw.py
+++ b/TurTLE/test/test_fftw.py
@@ -12,6 +12,11 @@ try:
 except:
     plt = None
 
+try:
+    import cmocean
+except:
+    cmocean = None
+
 def main():
     niterations = 10
     nlist = [16, 32, 48, 24, 64, 12]
@@ -32,10 +37,23 @@ def main():
                  sys.argv[1:])
         df = h5py.File(c.simname + '.h5', 'r')
         df = h5py.File(c.simname + '_fields.h5', 'r')
+        field0_complex = df['field0/complex/0'][...]
         field1_complex = df['field1/complex/0'][...]
         field1_real = df['field1/real/0'][...]
         npoints = field1_real.size//3
 
+        # first test symmetrize
+        # field0 is subject to TurTLE symmetrize
+        # field1 has additionally been through ift+dft+normalize
+        # ideally they should be identical. otherwise it means
+        # the TurTLE symmetrize is broken.
+        f0 = field0_complex[:, :, 0, :]
+        f1 = field1_complex[:, :, 0, :]
+        diff = np.abs(f0 - f1)
+        ii = np.unravel_index(np.argmax(diff), diff.shape)
+        assert(np.max(np.abs(diff)) < 1e-5)
+
+        # now compare numpy with fftw
         np_field1_real = np.fft.irfftn(field1_complex, axes = (0, 1, 2)).transpose(1, 0, 2, 3)
         L2normr = np.sqrt(np.mean(np.sum(field1_real**2, axis = 3)))
         np_L2normr = np.sqrt(np.mean(np.sum(np_field1_real**2, axis = 3)))
@@ -54,10 +72,24 @@ def main():
 
         if not type(plt) == type(None):
             f = plt.figure()
-            a = f.add_subplot(121)
+            a = f.add_subplot(221)
             a.imshow(np.log(np.abs(np_field1_complex[:, :, 0, 0])), interpolation = 'nearest')
-            a = f.add_subplot(122)
+            a.set_title('numpy')
+            a = f.add_subplot(222)
             a.imshow(np.log(np.abs(field1_complex[:, :, 0, 0])), interpolation = 'nearest')
+            a.set_title('TurTLE')
+            if not type(cmocean)  == type(None):
+                a = f.add_subplot(223)
+                a.imshow(
+                        np.log(np.angle(np_field1_complex[:, :, 0, 0])),
+                        interpolation = 'nearest',
+                        cmap  = cmocean.cm.phase)
+                a = f.add_subplot(224)
+                a.imshow(
+                        np.log(np.angle(field1_complex[:, :, 0, 0])),
+                        interpolation = 'nearest',
+                        cmap  = cmocean.cm.phase)
+            f.tight_layout()
             f.savefig(c.simname + '_complex_slice_kx0.pdf')
     return None
 
diff --git a/cpp/field.cpp b/cpp/field.cpp
index b92fce63..b647c98e 100644
--- a/cpp/field.cpp
+++ b/cpp/field.cpp
@@ -1455,9 +1455,9 @@ void field<rnumber, be, fc>::Hermitian_reflect()
 template <typename rnumber,
           field_backend be,
           field_components fc>
-void field<rnumber, be, fc>::symmetrize_FFTW()
+void field<rnumber, be, fc>::symmetrize_FFT()
 {
-    TIMEZONE("field::symmetrize_FFTW");
+    TIMEZONE("field::symmetrize_FFT");
     assert(!this->real_space_representation);
     typename fftw_interface<rnumber>::complex *cdata = this->get_cdata();
 
@@ -1501,7 +1501,6 @@ void field<rnumber, be, fc>::symmetrize_FFTW()
         }
     }
 
-    // for debugging, just use FFTW
     this->ift();
     this->dft();
     this->normalize();
@@ -1514,6 +1513,8 @@ template <typename rnumber,
 void field<rnumber, be, fc>::symmetrize()
 {
     TIMEZONE("field::symmetrize");
+    this->symmetrize_FFT();
+    return;
     assert(!this->real_space_representation);
     typename fftw_interface<rnumber>::complex *cdata = this->get_cdata();
 
@@ -1637,11 +1638,11 @@ void field<rnumber, be, fc>::symmetrize()
             for (ptrdiff_t iz = 1; iz < ptrdiff_t(this->clayout->sizes[1]); iz++)
             {
                 ptrdiff_t izz = (this->clayout->sizes[1] - iz);
-                ptrdiff_t cindex = this->get_cindex(0, iyy, izz);
+                ptrdiff_t cindex = this->get_cindex(0, iyy, iz);
                 for (int cc = 0; cc < int(ncomp(fc)); cc++)
                 {
-                    (*(cdata + ncomp(fc)*cindex + cc))[0] =  ((*(bufferp + ncomp(fc)*iz+cc))[0] + (*(bufferm + ncomp(fc)*iz+cc))[0])/2;
-                    (*(cdata + ncomp(fc)*cindex + cc))[1] =  ((*(bufferp + ncomp(fc)*iz+cc))[1] - (*(bufferm + ncomp(fc)*iz+cc))[1])/2;
+                    (*(cdata + ncomp(fc)*cindex + cc))[0] =  ((*(bufferp + ncomp(fc)*iz+cc))[0] + (*(bufferm + ncomp(fc)*izz+cc))[0])/2;
+                    (*(cdata + ncomp(fc)*cindex + cc))[1] =  ((*(bufferp + ncomp(fc)*iz+cc))[1] - (*(bufferm + ncomp(fc)*izz+cc))[1])/2;
                 }
             }
             ptrdiff_t cindex = this->get_cindex(0, iyy, 0);
@@ -1661,8 +1662,8 @@ void field<rnumber, be, fc>::symmetrize()
                 ptrdiff_t cindex = this->get_cindex(0, iyy, izz);
                 for (int cc = 0; cc < int(ncomp(fc)); cc++)
                 {
-                    (*(cdata + ncomp(fc)*cindex + cc))[0] =  ((*(bufferp + ncomp(fc)*iz+cc))[0] + (*(bufferm + ncomp(fc)*iz+cc))[0])/2;
-                    (*(cdata + ncomp(fc)*cindex + cc))[1] = -((*(bufferp + ncomp(fc)*iz+cc))[1] - (*(bufferm + ncomp(fc)*iz+cc))[1])/2;
+                    (*(cdata + ncomp(fc)*cindex + cc))[0] =  ((*(bufferp + ncomp(fc)*iz+cc))[0] + (*(bufferm + ncomp(fc)*izz+cc))[0])/2;
+                    (*(cdata + ncomp(fc)*cindex + cc))[1] = -((*(bufferp + ncomp(fc)*iz+cc))[1] - (*(bufferm + ncomp(fc)*izz+cc))[1])/2;
                 }
             }
             ptrdiff_t cindex = this->get_cindex(0, iyy, 0);
diff --git a/cpp/field.hpp b/cpp/field.hpp
index f9c88bcd..b8775f21 100644
--- a/cpp/field.hpp
+++ b/cpp/field.hpp
@@ -122,7 +122,7 @@ class field
         void ift();
         void normalize();
         void symmetrize();
-        void symmetrize_FFTW();
+        void symmetrize_FFT();
         void Hermitian_reflect();
 
         /* stats */
diff --git a/cpp/full_code/symmetrize_test.cpp b/cpp/full_code/symmetrize_test.cpp
index 111d3a83..b3058654 100644
--- a/cpp/full_code/symmetrize_test.cpp
+++ b/cpp/full_code/symmetrize_test.cpp
@@ -25,7 +25,6 @@
 
 #include <string>
 #include <cmath>
-#include <random>
 #include "symmetrize_test.hpp"
 #include "fftw_tools.hpp"
 #include "scope_timer.hpp"
@@ -62,6 +61,16 @@ int symmetrize_test<rnumber>::read_parameters()
     return EXIT_SUCCESS;
 }
 
+/** \brief Sanity check for symmetrize
+ *
+ * Generate a random field, apply symmetrize to it, call the result f0.
+ * Make a copy, and perform an ift+dft+normalize, call the result f1.
+ * Compare f0 with f1, they should be fairly close. Since we know that
+ * ift+dft+normalize only changes the field within numerical precision
+ * if it's properly Hermitian, then any difference between f0 and f1 is
+ * due to our own symmetrize being broken somehow.
+ * */
+
 template <typename rnumber>
 int symmetrize_test<rnumber>::do_work(void)
 {
@@ -79,9 +88,6 @@ int symmetrize_test<rnumber>::do_work(void)
             this->comm,
             fftw_planner_string_to_flag[this->fftw_plan_rigor]);
     DEBUG_MSG("finished allocating field1\n");
-    std::default_random_engine rgen;
-    std::normal_distribution<rnumber> rdist;
-    rgen.seed(1);
     kspace<FFTW,SMOOTH> *kk = new kspace<FFTW, SMOOTH>(
             test_field0->clayout, this->dkx, this->dky, this->dkz);
 
@@ -96,41 +102,14 @@ int symmetrize_test<rnumber>::do_work(void)
     }
 
     // fill up test_field0
-    *test_field0 = 0.0;
-    test_field0->real_space_representation = false;
-    kk->CLOOP_K2(
-            [&](ptrdiff_t cindex,
-                ptrdiff_t xindex,
-                ptrdiff_t yindex,
-                ptrdiff_t zindex,
-                double k2){
-            test_field0->cval(cindex, 0, 0) = rdist(rgen);
-            test_field0->cval(cindex, 0, 1) = rdist(rgen);
-            test_field0->cval(cindex, 1, 0) = rdist(rgen);
-            test_field0->cval(cindex, 1, 1) = rdist(rgen);
-            test_field0->cval(cindex, 2, 0) = rdist(rgen);
-            test_field0->cval(cindex, 2, 1) = rdist(rgen);
-            if (k2 > 0)
-            {
-                test_field0->cval(cindex, 0, 0) /= sqrt(k2);
-                test_field0->cval(cindex, 0, 1) /= sqrt(k2);
-                test_field0->cval(cindex, 1, 0) /= sqrt(k2);
-                test_field0->cval(cindex, 1, 1) /= sqrt(k2);
-                test_field0->cval(cindex, 2, 0) /= sqrt(k2);
-                test_field0->cval(cindex, 2, 1) /= sqrt(k2);
-            }
-            else
-            {
-                test_field0->cval(cindex, 0, 0) = 0;
-                test_field0->cval(cindex, 0, 1) = 0;
-                test_field0->cval(cindex, 1, 0) = 0;
-                test_field0->cval(cindex, 1, 1) = 0;
-                test_field0->cval(cindex, 2, 0) = 0;
-                test_field0->cval(cindex, 2, 1) = 0;
-            }
-            });
-    // dealias (?!)
-    kk->template dealias<rnumber, THREE>(test_field0->get_cdata());
+    *test_field0 = rnumber(0.0);
+    make_gaussian_random_field<rnumber, FFTW, THREE, SMOOTH>(
+            kk,
+            test_field0,
+            1,
+            -5./3,
+            kk->kM,
+            1.0);
     // make the field divergence free
     kk->template force_divfree<rnumber>(test_field0->get_cdata());
     // apply symmetrize to test_field0
@@ -140,13 +119,10 @@ int symmetrize_test<rnumber>::do_work(void)
     // make copy in test_field1
     // this MUST be made after symmetrizing test_field0
     // (alternatively, we may symmetrize test_field1 as well before the ift-dft cycle
-    test_field1->real_space_representation = false;
     *test_field1 = test_field0->get_cdata();
 
     // go back and forth with test_field1, to enforce symmetry
-    test_field1->ift();
-    test_field1->dft();
-    test_field1->normalize();
+    test_field1->symmetrize_FFT();
 
     // now compare the two fields
     double max_diff = 0;
@@ -196,6 +172,11 @@ int symmetrize_test<rnumber>::do_work(void)
     DEBUG_MSG("found maximum relative difference %g at ix = %ld, iy = %ld, iz = %ld, wavenumber = %g, amplitudes %g %g\n",
             max_diff, ix, iy, iz, k_at_max_diff, a0, a1);
 
+    test_field0->io(
+            this->simname + "_fields.h5",
+            "field0",
+            0,
+            false);
     test_field1->io(
             this->simname + "_fields.h5",
             "field1",
-- 
GitLab