Properly include negative frequencies due to convention change in AMReX

Author: MSABuschmann

source/gravitational_waves.cpp

@@ -68,6 +68,12 @@ void GravitationalWaves::ComputeSpectrum(
     modifier->SelectComponents(comps);
 
     for (int i = 0; i < 6; ++i) {
+#ifdef OLD_FFT
+        const amrex::BoxArray &ba = ld.boxArray();
+        const amrex::DistributionMapping &dm = ld.DistributionMap();
+        du_real[i].define(ba, dm, 1, 0);
+        du_imag[i].define(ba, dm, 1, 0);
+#endif
         utils::Fft(ld, comps[i] + idx_offset, du_real[i], du_imag[i],
                    sim->geom[lev], false, zero_padding);
     }
@@ -127,15 +133,31 @@ void GravitationalWaves::ComputeSpectrum(
             for (int b = lo.y; b <= hi.y; ++b) {
                 AMREX_PRAGMA_SIMD
                 for (int a = lo.x; a <= hi.x; ++a) {
+                    // To account of negative frequencies
+                    double multpl = (a == 0 || a == dimN / 2) ? 1. : 2.;
                     int abc[3] = {a, b, c};
                     double running_sum = 0;
 
+                    int li = a >= dimN / 2 ? a - dimN : a;
+                    int lj = b >= dimN / 2 ? b - dimN : b;
+                    int lk = c >= dimN / 2 ? c - dimN : c;
+                    unsigned int sq = li * li + lj * lj + lk * lk;
+                    unsigned long index;
+                    if (unbinned) {
+                        index = std::lower_bound(ks.begin(), ks.end(), sq) -
+                                ks.begin() + omp_get_thread_num() * kmax;
+                    } else {
+                        index = static_cast<long>(
+                            std::sqrt(static_cast<double>(sq)) + .5);
+                    }
+
                     for (int i = 0; i < 3; ++i) {
                         for (int j = 0; j < 3; ++j) {
                             for (int l = 0; l < 3; ++l) {
                                 for (int m = 0; m < 3; ++m) {
                                     running_sum +=
                                         GetLambda(i, j, l, m, abc, dimN) *
+                                        multpl *
                                         (du_real_arr[mat[i][j]]->operator()(
                                              a, b, c) *
                                              du_real_arr[mat[l][m]]->operator()(
@@ -148,20 +170,6 @@ void GravitationalWaves::ComputeSpectrum(
                             }
                         }
                     }
-
-                    int li = a >= dimN / 2 ? a - dimN : a;
-                    int lj = b >= dimN / 2 ? b - dimN : b;
-                    int lk = c >= dimN / 2 ? c - dimN : c;
-                    unsigned int sq = li * li + lj * lj + lk * lk;
-                    unsigned long index;
-                    if (unbinned) {
-                        index = std::lower_bound(ks.begin(), ks.end(), sq) -
-                                ks.begin() + omp_get_thread_num() * kmax;
-                    } else {
-                        index = static_cast<long>(
-                            std::sqrt(static_cast<double>(sq)) + .5);
-                    }
-
                     gw_spectrum[index] += running_sum;
                 }
             }

source/output_types/spectrum.cpp

@@ -73,14 +73,17 @@ void Spectrum::Compute(const int id, const hid_t file_id, Sledgehamr *sim) {
             for (int j = jl; j <= jh; ++j) {
                 AMREX_PRAGMA_SIMD
                 for (int k = kl; k <= kh; ++k) {
+                    // To account for negative frequencies
+                    double multpl = (i == 0 || i == dimN / 2) ? 1. : 2.;
                     int li = i >= dimN / 2 ? i - dimN : i;
                     int lj = j >= dimN / 2 ? j - dimN : j;
                     int lk = k >= dimN / 2 ? k - dimN : k;
                     unsigned int sq = li * li + lj * lj + lk * lk;
                     unsigned long index =
                         std::lower_bound(ks.begin(), ks.end(), sq) -
                         ks.begin() + omp_get_thread_num() * kmax;
-                    spectrum[index] += pre * field_fft_arr(i, j, k, 0) *
+                    spectrum[index] += multpl * pre *
+                                       field_fft_arr(i, j, k, 0) *
                                        field_fft_arr(i, j, k, 0);
                 }
             }