[broken] temporary state

cpp/full_code/bandpass_stats.cpp

@@ -54,8 +54,20 @@ int bandpass_stats<rnumber>::initialize(void)
     this->iteration_list = hdf5_tools::read_vector<int>(
             parameter_file,
             "/bandpass_stats/parameters/iteration_list");
-    this->max_velocity_estimate = hdf5_tools::read_value<double>(parameter_file, "/bandpass_stats/parameters/max_velocity_estimate");
+    this->max_velocity_estimate = hdf5_tools::read_value<double>(
+            parameter_file,
+            "/bandpass_stats/parameters/max_velocity_estimate");
+    this->k0list = hdf5_tools::read_vector<double>(
+            parameter_file,
+            "/bandpass_stats/parameters/k0list");
+    this->k1list = hdf5_tools::read_vector<double>(
+            parameter_file,
+            "/bandpass_stats/parameters/k1list");
+    this->filter_type = hdf5_tools::read_string(
+            parameter_file,
+            "/bandpass_stats/parameters/filter_type");
     H5Fclose(parameter_file);
+    assert(this->k0list.size() == this->k1list.size());
     if (this->myrank == 0)
     {
         // set caching parameters
@@ -100,18 +112,29 @@ int bandpass_stats<rnumber>::work_on_current_iteration(void)
             this->vorticity->fftw_plan_rigor);
 
 
+    /// compute the velocity, store in vel
     invert_curl(
         this->kk,
         this->vorticity,
         vel);
 
-    vel->ift();
+    /// allocate temporary filtered field container
+    field<rnumber, FFTW, THREE> *filtered_vel0 = new field<rnumber, FFTW, THREE>(
+            this->nx, this->ny, this->nz,
+            this->comm,
+            this->vorticity->fftw_plan_rigor);
+    field<rnumber, FFTW, THREE> *filtered_vel1 = new field<rnumber, FFTW, THREE>(
+            this->nx, this->ny, this->nz,
+            this->comm,
+            this->vorticity->fftw_plan_rigor);
 
+    /// make max vel estimate std::vector, required by compute_rspace_stats
     std::vector<double> max_vel_estimate;
     max_vel_estimate.resize(4, max_velocity_estimate / sqrt(3));
     max_vel_estimate[3] = max_velocity_estimate;
 
     /// initialize `stat_group`.
+    /// i.e. open HDF5 file for writing stats
     hid_t stat_group;
     if (this->myrank == 0)
         stat_group = H5Gopen(
@@ -121,16 +144,32 @@ int bandpass_stats<rnumber>::work_on_current_iteration(void)
     else
         stat_group = 0;
 
-    DEBUG_MSG("stat_group = %ld\n");
-    vel->compute_rspace_stats(
-            stat_group,
-            "velocity",
-            this->iteration / this->niter_out,
-            max_vel_estimate);
+    /// loop over all bands
+    for (int band_counter = 0; band_counter < this->k0list.size(); band_counter++)
+    {
+        *filtered_vel0 = *vel;
+        filtered_vel1->operator=(*vel);
+        filtered_vel1 =*vel;
+        this->kk->template filter_calibrated_ell<rnumber, THREE>(
+                this->filtered_vel0->get_cdata(),
+                4*acos(0) / this->k0list[band_counter],
+                this->filter_type);
+        this->kk->template filter_calibrated_ell<rnumber, THREE>(
+                this->filtered_vel1->get_cdata(),
+                4*acos(0) / this->k1list[band_counter],
+                this->filter_type);
+        filtered_vel->compute_rspace_stats(
+                stat_group,
+                "velocity_band" + std::to_string(band_counter),
+                this->iteration / this->niter_out,
+                max_vel_estimate);
+    }
 
     if (this->myrank == 0)
         H5Gclose(stat_group);
 
+    delete filtered_vel1;
+    delete filtered_vel0;
     delete vel;
     return EXIT_SUCCESS;
 }

cpp/full_code/bandpass_stats.hpp

@@ -44,6 +44,8 @@ class bandpass_stats: public NSVE_field_stats<rnumber>
         int niter_out;
         double max_velocity_estimate;
         kspace<FFTW, SMOOTH> *kk;
+        std::vector<double> k0list, k1list;
+        std::string filter_type;
 
         bandpass_stats(
                 const MPI_Comm COMMUNICATOR,