compute dealias factor on the fly

this approach implicitly solves the anisotropic grid problem

bfps/cpp/kspace.cpp

@@ -69,7 +69,7 @@ kspace<be, dt>::kspace(
                 if (ii <= int(this->layout->sizes[0]/2))
                     this->ky[i] = this->dky*ii;
                 else
-                    this->ky[i] = this->dky*(ii - int(this->layout->sizes[1]));
+                    this->ky[i] = this->dky*(ii - int(this->layout->sizes[0]));
             }
             for (i = 0; i<int(this->layout->sizes[1]); i++)
             {
@@ -116,8 +116,6 @@ kspace<be, dt>::kspace(
         std::fill_n(nshell_local, this->nshells, 0);
     });
 
-    std::vector<std::unordered_map<int, double>> dealias_filter_threaded(omp_get_max_threads());
-
     this->CLOOP_K2_NXMODES(
             [&](ptrdiff_t cindex,
                 ptrdiff_t xindex,
@@ -131,9 +129,6 @@ kspace<be, dt>::kspace(
                 kshell_local_thread.getMine()[int(knorm/this->dk)] += nxmodes*knorm;
                 nshell_local_thread.getMine()[int(knorm/this->dk)] += nxmodes;
             }
-            if (dt == SMOOTH){
-                dealias_filter_threaded[omp_get_thread_num()][int(round(k2 / this->dk2))] = exp(-36.0 * pow(k2/this->kM2, 18.));
-            }
     });
 
     // Merge results
@@ -141,14 +136,6 @@ kspace<be, dt>::kspace(
     kshell_local_thread.mergeParallel();
     nshell_local_thread.mergeParallel();
 
-    if (dt == SMOOTH){
-        for(int idxMerge = 0 ; idxMerge < int(dealias_filter_threaded.size()) ; ++idxMerge){
-            for(const auto kv : dealias_filter_threaded[idxMerge]){
-                this->dealias_filter[kv.first] = kv.second;
-            }
-        }
-    }
-
     MPI_Allreduce(
             nshell_local_thread.getMasterData(),
             &this->nshell.front(),
@@ -437,6 +424,7 @@ int kspace<be, dt>::filter_calibrated_ell(
         const double ell,
         std::string filter_type)
 {
+    TIMEZONE("kspace::filter_calibrated_ell");
     if (filter_type == std::string("sharp_Fourier_sphere"))
     {
         this->template low_pass<rnumber, fc>(
@@ -464,23 +452,22 @@ template <typename rnumber,
           field_components fc>
 void kspace<be, dt>::dealias(typename fftw_interface<rnumber>::complex *__restrict__ a)
 {
+    TIMEZONE("kspace::dealias");
     switch(dt)
     {
         case TWO_THIRDS:
             this->low_pass<rnumber, fc>(a, this->kM);
             break;
         case SMOOTH:
-            this->CLOOP_K2(
+            this->CLOOP(
                 [&](ptrdiff_t cindex,
                     ptrdiff_t xindex,
                     ptrdiff_t yindex,
-                    ptrdiff_t zindex,
-                    double k2){
-                    //double tval = this->dealias_filter[int(round(k2 / this->dk2))];
+                    ptrdiff_t zindex){
                     double kk2 = (pow(this->kx[xindex]/this->kMx, 2) +
                                   pow(this->ky[yindex]/this->kMy, 2) +
                                   pow(this->kz[zindex]/this->kMz, 2));
-                    double tval = exp(-36.0 * (pow(kk2, 18.)));
+                    double tval = exp(-36.0 * (pow(kk2, 18)));
                     for (unsigned int tcounter=0; tcounter<2*ncomp(fc); tcounter++)
                         ((rnumber*)a)[2*ncomp(fc)*cindex + tcounter] *= tval;
                 });

bfps/cpp/kspace.hpp

@@ -54,7 +54,6 @@ class kspace
         /* mode and dealiasing information */
         double kMx, kMy, kMz, kM, kM2;
         std::vector<double> kx, ky, kz;
-        std::unordered_map<int, double> dealias_filter;
         std::vector<double> kshell;
         std::vector<int64_t> nshell;
         int nshells;
@@ -156,8 +155,8 @@ class kspace
                         for (hsize_t xindex = 0; xindex < this->layout->subsizes[2]; xindex++)
                         {
                             double k2 = (this->kx[xindex]*this->kx[xindex] +
-                                  this->ky[yindex]*this->ky[yindex] +
-                                  this->kz[zindex]*this->kz[zindex]);
+                                         this->ky[yindex]*this->ky[yindex] +
+                                         this->kz[zindex]*this->kz[zindex]);
                             expression(cindex, xindex, yindex, zindex, k2);
                             cindex++;
                         }
@@ -179,7 +178,6 @@ class kspace
                                             + zindex*this->layout->subsizes[2];
                         hsize_t xindex = 0;
                         double k2 = (
-                                this->kx[xindex]*this->kx[xindex] +
                                 this->ky[yindex]*this->ky[yindex] +
                                 this->kz[zindex]*this->kz[zindex]);
                         expression(cindex, xindex, yindex, zindex, k2, 1);