diff --git a/cpp/full_code/bandpass_stats.cpp b/cpp/full_code/bandpass_stats.cpp index 938d4c0146b3a79625c9ed6119f01f5053d6bd7b..5154cfc2db1b469ee3345cad71abcf9a7214e1d3 100644 --- a/cpp/full_code/bandpass_stats.cpp +++ b/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; } diff --git a/cpp/full_code/bandpass_stats.hpp b/cpp/full_code/bandpass_stats.hpp index 05fae8e1a9ce5a900f44c09be02367d89d74e480..2158faf1a2233607bc4c539936439dfbd32797cb 100644 --- a/cpp/full_code/bandpass_stats.hpp +++ b/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,