Merge branch 'feature/static_field' into develop

CMakeLists.txt

@@ -186,6 +186,7 @@ set(cpp_for_lib
     ${PROJECT_SOURCE_DIR}/cpp/full_code/code_base.cpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/direct_numerical_simulation.cpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/NSVE.cpp
+    ${PROJECT_SOURCE_DIR}/cpp/full_code/static_field.cpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/joint_acc_vel_stats.cpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/test.cpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/filter_test.cpp
@@ -227,6 +228,7 @@ set(hpp_for_lib
     ${PROJECT_SOURCE_DIR}/cpp/full_code/code_base.hpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/direct_numerical_simulation.hpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/NSVE.hpp
+    ${PROJECT_SOURCE_DIR}/cpp/full_code/static_field.hpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/joint_acc_vel_stats.hpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/test.hpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/filter_test.hpp

TurTLE/DNS.py

@@ -440,7 +440,7 @@ class DNS(_code):
         assert (self.parameters['niter_todo'] % self.parameters['niter_stat'] == 0)
         assert (self.parameters['niter_todo'] % self.parameters['niter_out']  == 0)
         assert (self.parameters['niter_out']  % self.parameters['niter_stat'] == 0)
-        if self.dns_type in ['NSVEparticles_no_output', 'NSVEcomplex_particles', 'NSVEparticles']:
+        if self.dns_type in ['NSVEparticles_no_output', 'NSVEcomplex_particles', 'NSVEparticles', 'static_field']:
             assert (self.parameters['niter_todo'] % self.parameters['niter_part'] == 0)
             assert (self.parameters['niter_out']  % self.parameters['niter_part'] == 0)
         _code.write_par(self, iter0 = iter0)
@@ -643,6 +643,10 @@ class DNS(_code):
                 'NSVEparticles_no_output',
                 help = 'plain Navier-Stokes vorticity formulation, with basic fluid tracers, checkpoints are NOT SAVED')
 
+        parser_static_field = subparsers.add_parser(
+                'static_field',
+                help = 'static field with basic fluid tracers')
+
         parser_NSVEp2 = subparsers.add_parser(
                 'NSVEparticles',
                 help = 'plain Navier-Stokes vorticity formulation, with basic fluid tracers')
@@ -655,7 +659,7 @@ class DNS(_code):
                 'NSVEp_extra_sampling',
                 help = 'plain Navier-Stokes vorticity formulation, with basic fluid tracers, that sample velocity gradient, as well as pressure and its derivatives.')
 
-        for parser in ['NSVEparticles_no_output', 'NSVEp2', 'NSVEp2p', 'NSVEp_extra']:
+        for parser in ['NSVEparticles_no_output', 'NSVEp2', 'NSVEp2p', 'NSVEp_extra', 'static_field']:
             eval('self.simulation_parser_arguments({0})'.format('parser_' + parser))
             eval('self.job_parser_arguments({0})'.format('parser_' + parser))
             eval('self.particle_parser_arguments({0})'.format('parser_' + parser))
@@ -698,7 +702,7 @@ class DNS(_code):
         self.dns_type = opt.DNS_class
         self.name = self.dns_type + '-' + self.fluid_precision + '-v' + TurTLE.__version__
         # merge parameters if needed
-        if self.dns_type in ['NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output', 'NSVEp_extra_sampling']:
+        if self.dns_type in ['NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output', 'NSVEp_extra_sampling', 'static_field']:
             for k in self.NSVEp_extra_parameters.keys():
                 self.parameters[k] = self.NSVEp_extra_parameters[k]
         if type(extra_parameters) != type(None):
@@ -763,7 +767,7 @@ class DNS(_code):
             # hardcoded FFTW complex representation size
             field_size = 3*(opt.nx+2)*opt.ny*opt.nz*self.fluid_dtype.itemsize
             checkpoint_size = field_size
-            if self.dns_type in ['NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output', 'NSVEp_extra_sampling']:
+            if self.dns_type in ['static_field', 'NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output', 'NSVEp_extra_sampling']:
                 rhs_size = self.parameters['tracers0_integration_steps']
                 if type(opt.tracers0_integration_steps) != type(None):
                     rhs_size = opt.tracers0_integration_steps
@@ -1039,7 +1043,7 @@ class DNS(_code):
                 f['vorticity/complex/{0}'.format(0)] = data
             f.close()
         # now take care of particles' initial condition
-        if self.dns_type in ['NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output', 'NSVEp_extra_sampling']:
+        if self.dns_type in ['static_field', 'NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output', 'NSVEp_extra_sampling']:
             self.generate_particle_data(opt = opt)
         return None
     def launch_jobs(

cpp/full_code/direct_numerical_simulation.hpp

@@ -62,6 +62,15 @@ class direct_numerical_simulation: public code_base
         int read_iteration(void);
         int write_iteration(void);
         int grow_file_datasets(void);
+
+        const inline std::string get_current_fname()
+        {
+            return (
+                    std::string(this->simname) +
+                    std::string("_checkpoint_") +
+                    std::to_string(this->checkpoint) +
+                    std::string(".h5"));
+        }
 };
 
 #endif//DIRECT_NUMERICAL_SIMULATION_HPP

cpp/full_code/static_field.cpp

+/******************************************************************************
+*                                                                             *
+*  Copyright 2019 Max Planck Institute for Dynamics and Self-Organization     *
+*                                                                             *
+*  This file is part of bfps.                                                 *
+*                                                                             *
+*  bfps is free software: you can redistribute it and/or modify               *
+*  it under the terms of the GNU General Public License as published          *
+*  by the Free Software Foundation, either version 3 of the License,          *
+*  or (at your option) any later version.                                     *
+*                                                                             *
+*  bfps is distributed in the hope that it will be useful,                    *
+*  but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+*  GNU General Public License for more details.                               *
+*                                                                             *
+*  You should have received a copy of the GNU General Public License          *
+*  along with bfps.  If not, see <http://www.gnu.org/licenses/>               *
+*                                                                             *
+* Contact: Cristian.Lalescu@ds.mpg.de                                         *
+*                                                                             *
+******************************************************************************/
+
+
+
+#include <string>
+#include <cmath>
+#include "static_field.hpp"
+#include "scope_timer.hpp"
+#include "fftw_tools.hpp"
+
+
+template <typename rnumber>
+int static_field<rnumber>::initialize(void)
+{
+    TIMEZONE("static_file::initialize");
+    this->read_iteration();
+    this->read_parameters();
+    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 + ".h5").c_str(),
+                H5F_ACC_RDWR,
+                fapl);
+    }
+    int data_file_problem;
+    if (this->myrank == 0)
+        data_file_problem = this->grow_file_datasets();
+    MPI_Bcast(&data_file_problem, 1, MPI_INT, 0, this->comm);
+    if (data_file_problem > 0)
+    {
+        std::cerr <<
+            data_file_problem <<
+            " problems growing file datasets.\ntrying to exit now." <<
+            std::endl;
+        return EXIT_FAILURE;
+    }
+
+    this->vorticity = new field<rnumber, FFTW, THREE>(
+            nx, ny, nz,
+            this->comm,
+	    fftw_planner_string_to_flag[this->fftw_plan_rigor]);
+    this->vorticity->real_space_representation = false;
+
+    this->velocity = new field<rnumber, FFTW, THREE>(
+            nx, ny, nz,
+            this->comm,
+	    fftw_planner_string_to_flag[this->fftw_plan_rigor]);
+    this->velocity->real_space_representation = false;
+    //read vorticity field
+    this->vorticity->io(
+	    this->get_current_fname(),
+	    "vorticity",
+	    this->iteration,
+	    true);
+    this->kk = new kspace<FFTW, SMOOTH>(
+            this->vorticity->clayout, this->dkx, this->dky, this->dkz);
+
+    // compute the velocity field and store
+    invert_curl(
+        this->kk,
+        this->vorticity,
+        velocity);
+    // transform velocity and vorticity fields to real space
+    this->velocity->ift();
+    this->vorticity->ift();
+
+    // initialize particles
+    this->ps = particles_system_builder(
+                this->velocity,              // (field object)
+                this->kk,                     // (kspace object, contains dkx, dky, dkz)
+                tracers0_integration_steps, // to check coherency between parameters and hdf input file (nb rhs)
+                (long long int)nparticles,  // to check coherency between parameters and hdf input file
+	            this->get_current_fname(),    // particles input filename
+                std::string("/tracers0/state/") + std::to_string(this->iteration), // dataset name for initial input
+                std::string("/tracers0/rhs/")  + std::to_string(this->iteration),  // dataset name for initial input
+                tracers0_neighbours,        // parameter (interpolation no neighbours)
+                tracers0_smoothness,        // parameter
+                this->comm,
+                this->iteration+1);
+    // initialize output objects
+    this->particles_output_writer_mpi = new particles_output_hdf5<
+        long long int, double, 3>(
+                MPI_COMM_WORLD,
+                "tracers0",
+                nparticles,
+                tracers0_integration_steps);
+    this->particles_output_writer_mpi->setParticleFileLayout(this->ps->getParticleFileLayout());
+    this->particles_sample_writer_mpi = new particles_output_sampling_hdf5<
+        long long int, double, 3>(
+                MPI_COMM_WORLD,
+                this->ps->getGlobalNbParticles(),
+                (this->simname + "_particles.h5"),
+                "tracers0",
+                "position/0");
+    this->particles_sample_writer_mpi->setParticleFileLayout(this->ps->getParticleFileLayout());
+
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int static_field<rnumber>::step(void)
+{
+    TIMEZONE("static_file::step");
+    this->ps->completeLoop(this->dt);
+    this->iteration++;
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int static_field<rnumber>::write_checkpoint(void)
+{
+    TIMEZONE("static_file::write_checkpoint");
+    this->particles_output_writer_mpi->open_file(this->get_current_fname());
+    this->particles_output_writer_mpi->template save<3>(
+            this->ps->getParticlesState(),
+            this->ps->getParticlesRhs(),
+            this->ps->getParticlesIndexes(),
+            this->ps->getLocalNbParticles(),
+            this->iteration);
+    this->particles_output_writer_mpi->close_file();
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int static_field<rnumber>::finalize(void)
+{
+    TIMEZONE("static_field::finalize");
+    if (this->myrank == 0)
+        H5Fclose(this->stat_file);
+    // good practice rule number n+1: always delete in inverse order of allocation
+    delete this->ps.release();
+    delete this->particles_output_writer_mpi;
+    delete this->particles_sample_writer_mpi;
+    delete this->kk;
+    delete this->velocity;
+    delete this->vorticity;
+    return EXIT_SUCCESS;
+}
+
+/** \brief Compute statistics.
+ *
+ */
+
+template <typename rnumber>
+int static_field<rnumber>::do_stats()
+{
+    TIMEZONE("static_field::do_stats");
+    /// 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
+    if (!(this->iteration % this->niter_part == 0 ||
+          ((this->iteration + this->niter_part_fine_duration/2) % this->niter_part_fine_period <=
+           this->niter_part_fine_duration)))
+        return EXIT_SUCCESS;
+
+    // allocate temporary data array
+    std::unique_ptr<double[]> pdata(new double[3*this->ps->getLocalNbParticles()]);
+
+    /// copy position data
+
+    /// sample position
+    std::copy(this->ps->getParticlesState(),
+              this->ps->getParticlesState()+3*this->ps->getLocalNbParticles(),
+              pdata.get());
+    this->particles_sample_writer_mpi->template save_dataset<3>(
+            "tracers0",
+            "position",
+            this->ps->getParticlesState(),
+            &pdata,
+            this->ps->getParticlesIndexes(),
+            this->ps->getLocalNbParticles(),
+            this->ps->get_step_idx()-1);
+
+    /// sample velocity
+    std::fill_n(pdata.get(), 3*this->ps->getLocalNbParticles(), 0);
+    this->ps->sample_compute_field(*this->velocity, pdata.get());
+    this->particles_sample_writer_mpi->template save_dataset<3>(
+            "tracers0",
+            "velocity",
+            this->ps->getParticlesState(),
+            &pdata,
+            this->ps->getParticlesIndexes(),
+            this->ps->getLocalNbParticles(),
+            this->ps->get_step_idx()-1);
+
+    // deallocate temporary data array
+    delete[] pdata.release();
+
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int static_field<rnumber>::read_parameters(void)
+{
+    TIMEZONE("static_field::read_parameters");
+    this->direct_numerical_simulation::read_parameters();
+    hid_t parameter_file = H5Fopen((this->simname + ".h5").c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+    this->fftw_plan_rigor = hdf5_tools::read_string(parameter_file, "parameters/fftw_plan_rigor");
+    this->dt = hdf5_tools::read_value<double>(parameter_file, "parameters/dt");
+    this->niter_part = hdf5_tools::read_value<int>(parameter_file, "parameters/niter_part");
+    this->niter_part_fine_period = hdf5_tools::read_value<int>(parameter_file, "parameters/niter_part_fine_period");
+    this->niter_part_fine_duration = hdf5_tools::read_value<int>(parameter_file, "parameters/niter_part_fine_duration");
+    this->nparticles = hdf5_tools::read_value<int>(parameter_file, "parameters/nparticles");
+    this->tracers0_integration_steps = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_integration_steps");
+    this->tracers0_neighbours = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_neighbours");
+    this->tracers0_smoothness = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_smoothness");
+    H5Fclose(parameter_file);
+    return EXIT_SUCCESS;
+}
+
+template class static_field<float>;
+template class static_field<double>;
+

cpp/full_code/static_field.hpp

+/**********************************************************************
+*                                                                     *
+*  Copyright 2017 Max Planck Institute                                *
+*                 for Dynamics and Self-Organization                  *
+*                                                                     *
+*  This file is part of bfps.                                         *
+*                                                                     *
+*  bfps is free software: you can redistribute it and/or modify       *
+*  it under the terms of the GNU General Public License as published  *
+*  by the Free Software Foundation, either version 3 of the License,  *
+*  or (at your option) any later version.                             *
+*                                                                     *
+*  bfps is distributed in the hope that it will be useful,            *
+*  but WITHOUT ANY WARRANTY; without even the implied warranty of     *
+*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the      *
+*  GNU General Public License for more details.                       *
+*                                                                     *
+*  You should have received a copy of the GNU General Public License  *
+*  along with bfps.  If not, see <http://www.gnu.org/licenses/>       *
+*                                                                     *
+* Contact: Cristian.Lalescu@ds.mpg.de                                 *
+*                                                                     *
+**********************************************************************/
+
+
+
+#ifndef STATIC_FIELD_HPP
+#define STATIC_FIELD_HPP
+
+
+
+#include <cstdlib>
+#include "base.hpp"
+#include "kspace.hpp"
+#include "field.hpp"
+#include "full_code/direct_numerical_simulation.hpp"
+#include "particles/particles_system_builder.hpp"
+#include "particles/particles_output_hdf5.hpp"
+#include "particles/particles_sampling.hpp"
+
+template <typename rnumber>
+class static_field: public direct_numerical_simulation
+{
+    public:
+
+        /* parameters that are read in read_parameters */
+        double dt;
+        std::string fftw_plan_rigor;
+
+        /* particle parameters */
+        int niter_part;
+        int niter_part_fine_period;
+        int niter_part_fine_duration;
+        int nparticles;
+        int tracers0_integration_steps;
+        int tracers0_neighbours;
+        int tracers0_smoothness;
+
+        /* other stuff */
+        kspace<FFTW, SMOOTH> *kk;
+        field<rnumber, FFTW, THREE> *vorticity;
+        field<rnumber, FFTW, THREE> *velocity;
+
+        /* other stuff */
+        std::unique_ptr<abstract_particles_system<long long int, double>> ps;
+
+        particles_output_hdf5<long long int, double,3> *particles_output_writer_mpi;
+        particles_output_sampling_hdf5<long long int, double, 3> *particles_sample_writer_mpi;
+
+
+        static_field(
+                const MPI_Comm COMMUNICATOR,
+                const std::string &simulation_name):
+            direct_numerical_simulation(
+                    COMMUNICATOR,
+                    simulation_name){}
+        ~static_field(){}
+
+        int initialize(void);
+        int step(void);
+        int finalize(void);
+
+        virtual int read_parameters(void);
+        int write_checkpoint(void);
+        int do_stats(void);
+};
+
+#endif//STATIC_FIELD_HPP
+