diff --git a/TurTLE/test/test_fftw.py b/TurTLE/test/test_fftw.py index 1b1d4e6349a703d6119bcaa2db7354bb85e4b776..bb35935a4bad1660fd15893a162efd24ef7e3867 100644 --- a/TurTLE/test/test_fftw.py +++ b/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__': diff --git a/cpp/full_code/symmetrize_test.cpp b/cpp/full_code/symmetrize_test.cpp index b3058654f1e98067e7bf3bb82bd53114922de03e..c7f3d52618417012c7889bcdcf2f74b75bbef9ca 100644 --- a/cpp/full_code/symmetrize_test.cpp +++ b/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;