fix field io mechanism for Kraichnan model

TurTLE/DNS.py

@@ -680,6 +680,7 @@ class DNS(_code):
             dns_type):
         pars = {}
         if dns_type == 'kraichnan_field':
+            pars['output_velocity'] = int(1)
             pars['field_random_seed'] = int(1)
             pars['spectrum_slope'] = float(-5./3)
             pars['spectrum_k_cutoff'] = float(16)
@@ -1026,19 +1027,25 @@ class DNS(_code):
         # take care of fields' initial condition
         # first, check if initial field exists
         need_field = False
-        if not os.path.exists(self.get_checkpoint_0_fname()):
-            need_field = True
-        else:
-            f = h5py.File(self.get_checkpoint_0_fname(), 'r')
-            try:
-                dset = f['vorticity/complex/0']
-                need_field = (dset.shape == (self.parameters['ny'],
-                                             self.parameters['nz'],
-                                             self.parameters['nx']//2+1,
-                                             3))
-            except:
+        if self.dns_type in [
+                'static_field',
+                'NSVEparticles',
+                'NSVEcomplex_particles',
+                'NSVEparticles_no_output',
+                'NSVEp_extra_sampling']:
+            if not os.path.exists(self.get_checkpoint_0_fname()):
                 need_field = True
-            f.close()
+            else:
+                f = h5py.File(self.get_checkpoint_0_fname(), 'r')
+                try:
+                    dset = f['vorticity/complex/0']
+                    need_field = (dset.shape == (self.parameters['ny'],
+                                                 self.parameters['nz'],
+                                                 self.parameters['nx']//2+1,
+                                                 3))
+                except:
+                    need_field = True
+                f.close()
         if need_field:
             f = h5py.File(self.get_checkpoint_0_fname(), 'a')
             if len(opt.src_simname) > 0:
@@ -1065,6 +1072,11 @@ class DNS(_code):
                        amplitude = 0.05)
                 f['vorticity/complex/{0}'.format(0)] = data
             f.close()
+        if self.dns_type == 'kraichnan_field':
+            if not os.path.exists(self.get_checkpoint_0_fname()):
+                f = h5py.File(self.get_checkpoint_0_fname(), 'a')
+                f.create_group('velocity/real')
+                f.close()
         # now take care of particles' initial condition
         if self.dns_type in [
                 'kraichnan_field',
@@ -1081,6 +1093,9 @@ class DNS(_code):
         if not os.path.exists(self.get_data_file_name()):
             self.generate_initial_condition(opt = opt)
             self.write_par()
+        if (('test' in self.dns_type) or
+            (self.dns_type in ['kraichnan_field'])):
+            self.check_current_vorticity_exists = False
         self.run(
                 nb_processes = opt.nb_processes,
                 nb_threads_per_process = opt.nb_threads_per_process,

TurTLE/_code.py

@@ -176,6 +176,7 @@ class _code(_base):
                 """
         self.host_info = host_info
         self.main = ''
+        self.check_current_vorticity_exists = True
         return None
     def write_src(self):
         with open(self.name + '.cpp', 'w') as outfile:
@@ -278,7 +279,7 @@ class _code(_base):
         if os.path.exists(os.path.join(self.work_dir, 'stop_' + self.simname)):
             warnings.warn("Found stop_<simname> file, I will remove it.")
             os.remove(os.path.join(self.work_dir, 'stop_' + self.simname))
-        if not 'test' in self.name:
+        if self.check_current_vorticity_exists:
             with h5py.File(os.path.join(self.work_dir, self.simname + '.h5'), 'r') as data_file:
                 iter0 = data_file['iteration'][...]
                 checkpoint0 = data_file['checkpoint'][()]

cpp/field.cpp

@@ -2052,6 +2052,7 @@ int make_gaussian_random_field(
     rgen.resize(omp_get_max_threads());
     // seed random number generators such that no seed is ever repeated if we change the value of rseed.
     // basically use a multi-dimensional array indexing technique to come up with actual seed.
+    // Note: this method IS NOT MPI/OpenMP-invariant!
     for (int thread_id=0; thread_id < omp_get_max_threads(); thread_id++)
     {
         int current_seed = (

cpp/full_code/kraichnan_field.cpp

@@ -65,7 +65,7 @@ int kraichnan_field<rnumber>::initialize(void)
             nx, ny, nz,
             this->comm,
 	        fftw_planner_string_to_flag[this->fftw_plan_rigor]);
-    this->velocity->real_space_representation = false;
+    this->velocity->real_space_representation = true;
 
     // the velocity layout should be the same as the vorticity one
     this->kk = new kspace<FFTW, SMOOTH>(
@@ -104,30 +104,60 @@ int kraichnan_field<rnumber>::initialize(void)
     DEBUG_MSG("Coefficient: %g\n",
             this->spectrum_coefficient);
 
+    // is this the first iteration?
+    if ((this->iteration == 0) and (this->output_velocity == 1))
+    {
+        // if yes, generate random field and save it
+        this->generate_random_velocity();
+        this->velocity->io(
+                this->get_current_fname(),
+                "velocity",
+                this->iteration,
+                false);
+    }
+    else
+    {
+        // if not, read the random field that exists in the checkpoint file
+        this->velocity->io(
+                this->get_current_fname(),
+                "velocity",
+                this->iteration,
+                true);
+    }
+
     return EXIT_SUCCESS;
 }
 
 template <typename rnumber>
-int kraichnan_field<rnumber>::step(void)
+int kraichnan_field<rnumber>::generate_random_velocity(void)
 {
-    TIMEZONE("kraichnan_file::step");
-
+    TIMEZONE("kraichnan_file::make_random_field");
+    // generate Gaussian random field
     make_gaussian_random_field(
         this->kk,
         this->velocity,
-        this->iteration,
+        this->iteration, // not an ideal choice because resulting field sequence depends on MPI/OpenMP configuration
         this->spectrum_slope,
         this->spectrum_k_cutoff,
         this->spectrum_coefficient * 3./2.); // incompressibility projection factor
-    this->velocity->ift();
-
+    // this->velocity is now in Fourier space
+    // project divfree, requires field in Fourier space
     DEBUG_MSG("L2Norm before: %g\n",
             this->velocity->L2norm(this->kk));
     this->kk->template project_divfree<rnumber>(this->velocity->get_cdata());
     DEBUG_MSG("L2Norm after: %g\n",
             this->velocity->L2norm(this->kk));
+    // transform velocity to real space representation
+    this->velocity->ift();
+    return EXIT_SUCCESS;
+}
 
+template <typename rnumber>
+int kraichnan_field<rnumber>::step(void)
+{
+    TIMEZONE("kraichnan_file::step");
     this->ps->completeLoop(sqrt(this->dt));
+    this->generate_random_velocity();
     this->iteration++;
     return EXIT_SUCCESS;
 }
@@ -136,7 +166,19 @@ template <typename rnumber>
 int kraichnan_field<rnumber>::write_checkpoint(void)
 {
     TIMEZONE("kraichnan_file::write_checkpoint");
-    this->io_checkpoint(false);
+    this->update_checkpoint();
+    // output field
+    if (this->output_velocity == 1)
+    {
+        assert(this->velocity->real_space_representation);
+        std::string fname = this->get_current_fname();
+        this->velocity->io(
+                fname,
+                "velocity",
+                this->iteration,
+                false);
+    }
+    // output particles
     this->particles_output_writer_mpi->open_file(this->get_current_fname());
     this->particles_output_writer_mpi->template save<3>(
             this->ps->getParticlesState(),
@@ -172,6 +214,32 @@ template <typename rnumber>
 int kraichnan_field<rnumber>::do_stats()
 {
     TIMEZONE("kraichnan_field::do_stats");
+    /// compute field statistics
+    if (this->iteration % this->niter_stat == 0)
+    {
+        hid_t stat_group;
+        if (this->myrank == 0)
+            stat_group = H5Gopen(
+                    this->stat_file,
+                    "statistics",
+                    H5P_DEFAULT);
+        else
+            stat_group = 0;
+        field<rnumber, FFTW, THREE> *tvelocity = new field<rnumber, FFTW, THREE>(
+                this->nx, this->ny, this->nz,
+                this->comm,
+	            fftw_planner_string_to_flag[this->fftw_plan_rigor]);
+        *tvelocity = *this->velocity;
+        tvelocity->compute_stats(
+            this->kk,
+            stat_group,
+            "velocity",
+            this->iteration / this->niter_stat,
+            this->max_velocity_estimate/sqrt(3));
+        if (this->myrank == 0)
+            H5Gclose(stat_group);
+        delete tvelocity;
+    }
     /// either one of two conditions suffices to compute statistics:
     /// 1) current iteration is a multiple of niter_part
     /// 2) we are within niter_part_fine_duration/2 of a multiple of niter_part_fine_period
@@ -235,6 +303,8 @@ int kraichnan_field<rnumber>::read_parameters(void)
     this->spectrum_k_cutoff = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_k_cutoff");
     this->spectrum_coefficient = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_coefficient");
     this->field_random_seed = hdf5_tools::read_value<int>(parameter_file, "parameters/field_random_seed");
+    this->output_velocity = hdf5_tools::read_value<int>(parameter_file, "parameters/output_velocity");
+    this->max_velocity_estimate = hdf5_tools::read_value<double>(parameter_file, "parameters/max_velocity_estimate");
     H5Fclose(parameter_file);
     return EXIT_SUCCESS;
 }
@@ -242,6 +312,8 @@ int kraichnan_field<rnumber>::read_parameters(void)
 template <class rnumber>
 void kraichnan_field<rnumber>::update_checkpoint()
 {
+    TIMEZONE("kraichnan_field::update_checkpoint");
+    DEBUG_MSG("entered kraichan_field::update_checkpoint\n");
     std::string fname = this->get_current_fname();
     if (this->kk->layout->myrank == 0)
     {
@@ -297,6 +369,7 @@ void kraichnan_field<rnumber>::update_checkpoint()
         }
     }
     MPI_Bcast(&this->checkpoint, 1, MPI_INT, 0, this->kk->layout->comm);
+    DEBUG_MSG("exiting kraichan_field::update_checkpoint\n");
 }
 
 template class kraichnan_field<float>;

cpp/full_code/kraichnan_field.hpp

@@ -55,6 +55,7 @@ class kraichnan_field: public direct_numerical_simulation
         int tracers0_integration_steps;
         int tracers0_neighbours;
         int tracers0_smoothness;
+        int output_velocity;
 
         /* other stuff */
         kspace<FFTW, SMOOTH> *kk;
@@ -62,6 +63,7 @@ class kraichnan_field: public direct_numerical_simulation
         double spectrum_slope;
         double spectrum_k_cutoff;
         double spectrum_coefficient;
+        double max_velocity_estimate;
         int field_random_seed;
 
         /* other stuff */
@@ -87,20 +89,8 @@ class kraichnan_field: public direct_numerical_simulation
         int write_checkpoint(void);
         int do_stats(void);
 
-        
+        int generate_random_velocity(void);
         void update_checkpoint(void);
-        inline void io_checkpoint(bool read = true)
-        {
-            assert(this->velocity->real_space_representation);
-            if (!read)
-                this->update_checkpoint();
-            std::string fname = this->get_current_fname();
-            this->velocity->io(
-                    fname,
-                    "velocity",
-                    this->iteration,
-                    read);
-        }
 };
 
 #endif//KRAICHNAN_FIELD_HPP