[wip] adds code for Rij(x, 0, 0) etc extraction

TurTLE/PP.py

@@ -104,6 +104,8 @@ class PP(_code):
             self,
             dns_type = 'joint_acc_vel_stats'):
         pars = {}
+        if dns_type == 'get_3D_correlations':
+            pars['full_snapshot_output'] = int(False)
         if dns_type == 'pressure_stats':
             pars['max_pressure_estimate'] = float(1)
             pars['histogram_bins'] = int(129)

cpp/full_code/get_3D_correlations.cpp

@@ -71,9 +71,41 @@ int get_3D_correlations<rnumber>::initialize(void)
     this->iteration_list = hdf5_tools::read_vector<int>(
             parameter_file,
             "/get_3D_correlations/parameters/iteration_list");
+    this->full_snapshot_output = bool(
+            hdf5_tools::read_value<int>(
+                parameter_file,
+                "/get_3D_correlations/parameters/full_snapshot_output"));
     H5Fclose(parameter_file);
     // the following ensures the file is free for rank 0 to open in read/write mode
     MPI_Barrier(this->comm);
+    if (this->myrank == 0) {
+        // set caching parameters
+        hid_t fapl = H5Pcreate(H5P_FILE_ACCESS);
+        herr_t cache_err = H5Pset_cache(fapl, 0, 521, 134217728, 1.0);
+        variable_used_only_in_assert(cache_err);
+        DEBUG_MSG("when setting stat_file cache I got %d\n", cache_err);
+        this->stat_file = H5Fopen(
+                (this->simname + "_post.h5").c_str(),
+                H5F_ACC_RDWR,
+                fapl);
+    } else {
+        this->stat_file = 0;
+        }
+    int data_file_problem;
+    if (this->myrank == 0)
+        data_file_problem = hdf5_tools::require_size_file_datasets(
+                this->stat_file,
+                "get_3D_correlations",
+                (this->iteration_list.back() / this->niter_out) + 1);
+    MPI_Bcast(&data_file_problem, 1, MPI_INT, 0, this->comm);
+    if (data_file_problem > 0)
+    {
+        std::cerr <<
+            data_file_problem <<
+            " problems setting sizes of file datasets.\ntrying to exit now." <<
+            std::endl;
+        return EXIT_FAILURE;
+    }
     return EXIT_SUCCESS;
 }
 
@@ -81,12 +113,14 @@ template <typename rnumber>
 int compute_correlation(
         kspace<FFTW, SMOOTH> *kk,
         field<rnumber, FFTW, THREE> *vel,
-        field<rnumber, FFTW, THREExTHREE> *Rij)
+        field<rnumber, FFTW, THREExTHREE> *Rij,
+        hid_t group,
+        const std::string field_fname = "",
+        const int iteration = 0)
 {
     *(Rij) = 0.0;
     Rij->real_space_representation = false;
 
-
     const double dkvolume = kk->dkx*kk->dky*kk->dkz;
     // compute Rij
     kk->CLOOP_K2(
@@ -110,6 +144,49 @@ int compute_correlation(
 
     // go to real space
     Rij->ift();
+
+    /// output Rij field
+    if (field_fname.size() > 0) {
+        std::cerr << "fname size is " << field_fname.size()
+            << ", and fname is " << field_fname << std::endl;
+        Rij->io(
+                field_fname,
+                "Rij",
+                iteration,
+                false);
+    }
+
+    // extract angle-averaged data
+    double Rijx[9*Rij->get_nx()];
+    double Rijy[9*Rij->get_ny()];
+    double Rijz[9*Rij->get_nz()];
+    double Rijz_local[9*Rij->get_nz()];
+    std::fill_n(Rijz_local, 9*Rij->get_nz(), 0);
+    for (int i=0; i<3; i++)
+    {
+        for (int j=0; j<3; j++)
+        {
+            if (Rij->rlayout->myrank == Rij->rlayout->rank[0][0])
+            {
+                for (int ix=0; ix<Rij->get_nx(); ix++)
+                    Rijx[(i*3+j)*Rij->get_nx() + ix] = Rij->rval(
+                            ix, i, j);
+                for (int iy=0; iy<Rij->get_ny(); iy++)
+                    Rijy[(i*3+j)*Rij->get_ny() + iy] = Rij->rval(
+                            iy*(Rij->get_nx()+2), i, j);
+            }
+            for (int iiz=0; iiz<Rij->rlayout->subsizes[2]; iiz++)
+                Rijz_local[(i*3+j)*Rij->get_nz() + iiz + Rij->rlayout->starts[0]] =
+                    Rij->rval(iiz*(Rij->get_nx()+2)*Rij->get_ny(), i, j);
+        }
+    }
+    MPI_Allreduce(
+            Rijz_local,
+            Rijz,
+            9*Rij->get_nz(),
+            MPI_DOUBLE,
+            MPI_SUM,
+            Rij->comm);
     return EXIT_SUCCESS;
 }
 
@@ -124,20 +201,35 @@ int get_3D_correlations<rnumber>::work_on_current_iteration(void)
         this->vorticity,
         this->vel);
 
-    compute_correlation(this->kk, this->vel, this->Rij);
-
-    const std::string fname = (
+   std::string fname = "";
+   if (this->full_snapshot_output) {
+       fname = (
             this->simname +
             std::string("_checkpoint_Rij_") +
             std::to_string(this->iteration / (this->niter_out*this->checkpoints_per_file)) +
             std::string(".h5"));
+   }
+    /// initialize `stat_group`.
+    hid_t stat_group;
+    if (this->myrank == 0)
+        stat_group = H5Gopen(
+                this->stat_file,
+                "get_3D_correlations",
+                H5P_DEFAULT);
+    else
+        stat_group = 0;
 
-    /// output Rij field
-    this->Rij->io(
+    compute_correlation(
+            this->kk,
+            this->vel,
+            this->Rij,
+            stat_group,
             fname,
-            "Rij",
-            this->iteration,
-            false);
+            this->iteration);
+
+    /// close stat group
+    if (this->myrank == 0)
+        H5Gclose(stat_group);
 
     return EXIT_SUCCESS;
 }
@@ -149,6 +241,8 @@ int get_3D_correlations<rnumber>::finalize(void)
     delete this->Rij;
     delete this->vel;
     delete this->kk;
+    if (this->myrank == 0)
+        H5Fclose(this->stat_file);
     this->NSVE_field_stats<rnumber>::finalize();
     return EXIT_SUCCESS;
 }

cpp/full_code/get_3D_correlations.hpp

@@ -32,8 +32,10 @@ template <typename rnumber>
 class get_3D_correlations: public NSVE_field_stats<rnumber>
 {
     public:
+        hid_t stat_file;
         int checkpoints_per_file;
         int niter_out;
+        bool full_snapshot_output;
         kspace<FFTW, SMOOTH> *kk;
         field<rnumber, FFTW, THREE> *vel;
         field<rnumber, FFTW, THREExTHREE> *Rij;