adds hard symmetrize test

TurTLE/test/test_fftw.py

@@ -51,6 +51,7 @@ def main():
         f1 = field1_complex[:, :, 0, :]
         diff = np.abs(f0 - f1)
         ii = np.unravel_index(np.argmax(diff), diff.shape)
+        print('found maximum difference {0} between field0 and field1 at {1}'.format(diff[ii[0], ii[1], ii[2]], ii))
         assert(np.max(np.abs(diff)) < 1e-5)
 
         # now compare numpy with fftw
@@ -58,6 +59,7 @@ def main():
         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)))
         err = np.max(np.abs(field1_real - np_field1_real*npoints)) / L2normr
+        print('found maximum difference {0} between field1_real and np_field1_real'.format(err))
         assert(err < 1e-5)
 
         np_field1_complex = np.fft.rfftn(field1_real.transpose(1, 0, 2, 3), axes = (0, 1, 2)) / npoints
@@ -91,6 +93,20 @@ def main():
                         cmap  = cmocean.cm.phase)
             f.tight_layout()
             f.savefig(c.simname + '_complex_slice_kx0.pdf')
+
+        # hard test.
+        # generate random rseed2_field
+        # field0 is rseed2_field through TurTLE symmetrize
+        # field1 is rseed2_field through ift+dft+normalize
+        # compare field0 with field1
+        field0_complex = df['rseed2_field0/complex/0'][...]
+        field1_complex = df['rseed2_field1/complex/0'][...]
+        f0 = field0_complex[:, :, 0, :]
+        f1 = field1_complex[:, :, 0, :]
+        diff = np.abs(f0 - f1)
+        ii = np.unravel_index(np.argmax(diff), diff.shape)
+        print('found maximum difference {0} between rseed2_field0 and rseed2_field1 at {1}'.format(diff[ii[0], ii[1], ii[2]], ii))
+        assert(np.max(np.abs(diff)) < 1e-5)
     return None
 
 if __name__ == '__main__':

cpp/full_code/symmetrize_test.cpp

@@ -101,6 +101,10 @@ int symmetrize_test<rnumber>::do_work(void)
         H5Fclose(stat_file);
     }
 
+    /*************************************************
+     * here we test that after applying symmetrize
+     * field is actually Hermitian-symmetric
+     *************************************************/
     // fill up test_field0
     *test_field0 = rnumber(0.0);
     make_gaussian_random_field<rnumber, FFTW, THREE, SMOOTH>(
@@ -124,12 +128,15 @@ int symmetrize_test<rnumber>::do_work(void)
     // go back and forth with test_field1, to enforce symmetry
     test_field1->symmetrize_FFT();
 
-    // now compare the two fields
-    double max_diff = 0;
+    // temporary variables for field comparison
+    double max_diff;
     ptrdiff_t ix, iy, iz;
-    double k_at_max_diff = 0;
+    double k_at_max_diff;
     double a0, a1;
 
+    // now compare the two fields
+    max_diff = 0;
+    k_at_max_diff = 0;
     kk->CLOOP_K2(
             [&](ptrdiff_t cindex,
                 ptrdiff_t xindex,
@@ -169,7 +176,7 @@ int symmetrize_test<rnumber>::do_work(void)
                 a1 = sqrt(amplitude1);
             }
             });
-    DEBUG_MSG("found maximum relative difference %g at ix = %ld, iy = %ld, iz = %ld, wavenumber = %g, amplitudes %g %g\n",
+    DEBUG_MSG("rseed1 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(
@@ -189,6 +196,86 @@ int symmetrize_test<rnumber>::do_work(void)
             0,
             false);
 
+    /*************************************************
+     * here we test that applying our symmetrize, or
+     * just using fft, leads to the same result.
+     * there's no guarantee that this will work though,
+     * since FFT documentation doesn't say IFT+DFT+normalize leads to
+     * arithmetic mean-based Hermitian symmetry.
+     *************************************************/
+
+    make_gaussian_random_field<rnumber, FFTW, THREE, SMOOTH>(
+            kk,
+            test_field0,
+            2,
+            -5./3,
+            kk->kM,
+            1.0);
+
+    // copy unmodified data
+    *test_field1 = test_field0->get_cdata();
+
+    // apply different symmetrize methods
+    test_field0->symmetrize();
+    test_field1->symmetrize_FFT();
+
+    // now compare the two fields
+    max_diff = 0;
+    k_at_max_diff = 0;
+    kk->CLOOP_K2(
+            [&](ptrdiff_t cindex,
+                ptrdiff_t xindex,
+                ptrdiff_t yindex,
+                ptrdiff_t zindex,
+                double k2){
+            double diff_re0 = test_field0->cval(cindex, 0, 0) - test_field1->cval(cindex, 0, 0);
+            double diff_re1 = test_field0->cval(cindex, 1, 0) - test_field1->cval(cindex, 1, 0);
+            double diff_re2 = test_field0->cval(cindex, 2, 0) - test_field1->cval(cindex, 2, 0);
+            double diff_im0 = test_field0->cval(cindex, 0, 1) - test_field1->cval(cindex, 0, 1);
+            double diff_im1 = test_field0->cval(cindex, 1, 1) - test_field1->cval(cindex, 1, 1);
+            double diff_im2 = test_field0->cval(cindex, 2, 1) - test_field1->cval(cindex, 2, 1);
+            double diff = sqrt(diff_re0*diff_re0 + diff_re1*diff_re1 + diff_re2*diff_re2 +
+                               diff_im0*diff_im0 + diff_im1*diff_im1 + diff_im2*diff_im2);
+            double amplitude0 = (test_field0->cval(cindex, 0, 0)*test_field0->cval(cindex, 0, 0) +
+                                 test_field0->cval(cindex, 1, 0)*test_field0->cval(cindex, 1, 0) +
+                                 test_field0->cval(cindex, 2, 0)*test_field0->cval(cindex, 2, 0) +
+                                 test_field0->cval(cindex, 0, 1)*test_field0->cval(cindex, 0, 1) +
+                                 test_field0->cval(cindex, 1, 1)*test_field0->cval(cindex, 1, 1) +
+                                 test_field0->cval(cindex, 2, 1)*test_field0->cval(cindex, 2, 1));
+            double amplitude1 = (test_field1->cval(cindex, 0, 0)*test_field1->cval(cindex, 0, 0) +
+                                 test_field1->cval(cindex, 1, 0)*test_field1->cval(cindex, 1, 0) +
+                                 test_field1->cval(cindex, 2, 0)*test_field1->cval(cindex, 2, 0) +
+                                 test_field1->cval(cindex, 0, 1)*test_field1->cval(cindex, 0, 1) +
+                                 test_field1->cval(cindex, 1, 1)*test_field1->cval(cindex, 1, 1) +
+                                 test_field1->cval(cindex, 2, 1)*test_field1->cval(cindex, 2, 1));
+            double amplitude = sqrt((amplitude0 + amplitude1)/2);
+            if (amplitude > 0)
+            if (diff/amplitude > max_diff)
+            {
+                max_diff = diff / amplitude;
+                ix = xindex;
+                iy = yindex + test_field0->clayout->starts[0];
+                iz = zindex;
+                k_at_max_diff = sqrt(k2);
+                a0 = sqrt(amplitude0);
+                a1 = sqrt(amplitude1);
+            }
+            });
+    DEBUG_MSG("rseed2 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);
+
+    // output fields for possible subsequent comparison
+    test_field0->io(
+            this->simname + "_fields.h5",
+            "rseed2_field0",
+            0,
+            false);
+    test_field1->io(
+            this->simname + "_fields.h5",
+            "rseed2_field1",
+            0,
+            false);
+
     // deallocate
     delete kk;
     delete test_field1;