Merge branch 'feature/new-wrapper' into develop

Add simplified wrapper called `DNS`.
While this will make it slightly more complicated to generate
peculiar custom codes, the experience so far has been that this current
version is quite adequate for our needs.

Most of the code generation is removed from the Python wrapper, and
plain old C++ classes are written for each type of simulation that is to
be run.

bfps/__main__.py

@@ -28,6 +28,7 @@ import sys
 import argparse
 
 import bfps
+from .DNS import DNS
 from .NavierStokes import NavierStokes
 from .NSVorticityEquation import NSVorticityEquation
 from .FluidResize import FluidResize
@@ -64,13 +65,22 @@ def main():
                'NSManyParticles-double']
     parser.add_argument(
             'base_class',
-            choices = NSoptions + NSVEoptions + FRoptions + FCoptions + NSMPopt,
+            choices = ['DNS'] +
+                      NSoptions +
+                      NSVEoptions +
+                      FRoptions +
+                      FCoptions +
+                      NSMPopt,
             type = str)
     # first option is the choice of base class or -h or -v
     # all other options are passed on to the base_class instance
     opt = parser.parse_args(sys.argv[1:2])
     # error is thrown if first option is not a base class, so launch
     # cannot be executed by mistake.
+    if opt.base_class == 'DNS':
+        c = DNS()
+        c.launch(args = sys.argv[2:])
+        return None
     if 'double' in opt.base_class:
         precision = 'double'
     else:

bfps/_base.py

@@ -74,18 +74,28 @@ class _base(object):
             self,
             parameters = None,
             function_suffix = '',
-            file_group = 'parameters'):
+            template_class = '',
+            template_prefix = '',
+            file_group = 'parameters',
+            just_declaration = False,
+            simname_variable = 'simname'):
         if type(parameters) == type(None):
             parameters = self.parameters
         key = sorted(list(parameters.keys()))
-        src_txt = ('int read_parameters' + function_suffix + '()\n{\n' +
-                   'hid_t parameter_file;\n' +
-                   'hid_t dset, memtype, space;\n' +
-                   'char fname[256];\n' +
-                   'hsize_t dims[1];\n' +
-                   'char *string_data;\n' +
-                   'sprintf(fname, "%s.h5", simname);\n' +
-                   'parameter_file = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);\n')
+        src_txt = (template_prefix +
+                   'int ' +
+                   template_class +
+                   'read_parameters' + function_suffix + '()')
+        if just_declaration:
+            return src_txt + ';\n'
+        src_txt += ('\n{\n' +
+                    'hid_t parameter_file;\n' +
+                    'hid_t dset, memtype, space;\n' +
+                    'char fname[256];\n' +
+                    'hsize_t dims[1];\n' +
+                    'char *string_data;\n' +
+                    'sprintf(fname, "%s.h5", {0});\n'.format(simname_variable) +
+                    'parameter_file = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);\n')
         for i in range(len(key)):
             src_txt += 'dset = H5Dopen(parameter_file, "/{0}/{1}", H5P_DEFAULT);\n'.format(
                     file_group, key[i])

bfps/cpp/base.hpp

@@ -24,6 +24,7 @@
 
 
 
+#include <cassert>
 #include <mpi.h>
 #include <stdarg.h>
 #include <iostream>

bfps/cpp/full_code/NSVE.cpp

+#include <string>
+#include <cmath>
+#include "NSVE.hpp"
+
+
+template <typename rnumber>
+int NSVE<rnumber>::initialize(void)
+{
+    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);
+        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->fs = new vorticity_equation<rnumber, FFTW>(
+            simname.c_str(),
+            nx, ny, nz,
+            dkx, dky, dkz,
+            FFTW_MEASURE);
+    this->tmp_vec_field = new field<rnumber, FFTW, THREE>(
+            nx, ny, nz,
+            this->comm,
+            FFTW_MEASURE);
+
+
+    this->fs->checkpoints_per_file = checkpoints_per_file;
+    this->fs->nu = nu;
+    this->fs->fmode = fmode;
+    this->fs->famplitude = famplitude;
+    this->fs->fk0 = fk0;
+    this->fs->fk1 = fk1;
+    strncpy(this->fs->forcing_type, forcing_type, 128);
+    this->fs->iteration = this->iteration;
+    this->fs->checkpoint = this->checkpoint;
+
+    this->fs->cvorticity->real_space_representation = false;
+    this->fs->io_checkpoint();
+
+    if (this->myrank == 0 && this->iteration == 0)
+        this->fs->kk->store(stat_file);
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVE<rnumber>::main_loop(void)
+{
+    clock_t time0, time1;
+    double time_difference, local_time_difference;
+    time0 = clock();
+    bool stop_code_now = false;
+    int max_iter = (this->iteration + this->niter_todo -
+                    (this->iteration % this->niter_todo));
+    for (; this->iteration < max_iter;)
+    {
+    #ifdef USE_TIMINGOUTPUT
+        const std::string loopLabel = ("code::main_start::loop-" +
+                                       std::to_string(this->iteration));
+        TIMEZONE(loopLabel.c_str());
+    #endif
+        if (this->iteration % this->niter_stat == 0)
+            this->do_stats();
+        this->fs->step(dt);
+        this->iteration = this->fs->iteration;
+        if (this->fs->iteration % this->niter_out == 0)
+        {
+            this->fs->io_checkpoint(false);
+            this->checkpoint = this->fs->checkpoint;
+            this->write_iteration();
+        }
+        if (stop_code_now)
+            break;
+        time1 = clock();
+        local_time_difference = ((
+                (unsigned int)(time1 - time0)) /
+                ((double)CLOCKS_PER_SEC));
+        time_difference = 0.0;
+        MPI_Allreduce(
+                &local_time_difference,
+                &time_difference,
+                1,
+                MPI_DOUBLE,
+                MPI_SUM,
+                MPI_COMM_WORLD);
+        if (this->myrank == 0)
+            std::cout << "iteration " << iteration << " took " << time_difference/nprocs << " seconds" << std::endl;
+        if (this->myrank == 0)
+            std::cerr << "iteration " << iteration << " took " << time_difference/nprocs << " seconds" << std::endl;
+        time0 = time1;
+    }
+    if (this->iteration % this->niter_stat == 0)
+        this->do_stats();
+    time1 = clock();
+    local_time_difference = ((
+            (unsigned int)(time1 - time0)) /
+            ((double)CLOCKS_PER_SEC));
+    time_difference = 0.0;
+    MPI_Allreduce(
+            &local_time_difference,
+            &time_difference,
+            1,
+            MPI_DOUBLE,
+            MPI_SUM,
+            MPI_COMM_WORLD);
+    if (this->myrank == 0)
+        std::cout << "iteration " << iteration << " took " << time_difference/nprocs << " seconds" << std::endl;
+    if (this->myrank == 0)
+        std::cerr << "iteration " << iteration << " took " << time_difference/nprocs << " seconds" << std::endl;
+    if (this->iteration % this->niter_out != 0)
+    {
+        this->fs->io_checkpoint(false);
+        this->checkpoint = this->fs->checkpoint;
+        this->write_iteration();
+    }
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVE<rnumber>::finalize(void)
+{
+    if (this->myrank == 0)
+        H5Fclose(this->stat_file);
+    delete this->fs;
+    delete this->tmp_vec_field;
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVE<rnumber>::do_stats()
+{
+    hid_t stat_group;
+    if (this->myrank == 0)
+        stat_group = H5Gopen(
+                this->stat_file,
+                "statistics",
+                H5P_DEFAULT);
+    else
+        stat_group = 0;
+    fs->compute_velocity(fs->cvorticity);
+    *tmp_vec_field = fs->cvelocity->get_cdata();
+    tmp_vec_field->compute_stats(
+            fs->kk,
+            stat_group,
+            "velocity",
+            fs->iteration / niter_stat,
+            max_velocity_estimate/sqrt(3));
+
+    *tmp_vec_field = fs->cvorticity->get_cdata();
+    tmp_vec_field->compute_stats(
+            fs->kk,
+            stat_group,
+            "vorticity",
+            fs->iteration / niter_stat,
+            max_vorticity_estimate/sqrt(3));
+
+    if (this->myrank == 0)
+        H5Gclose(stat_group);
+
+    if (myrank == 0)
+    {
+        std::string fname = (
+                std::string("stop_") +
+                std::string(simname));
+        {
+            struct stat file_buffer;
+            this->stop_code_now = (
+                    stat(fname.c_str(), &file_buffer) == 0);
+        }
+    }
+    MPI_Bcast(
+            &this->stop_code_now,
+            1,
+            MPI_C_BOOL,
+            0,
+            MPI_COMM_WORLD);
+    return EXIT_SUCCESS;
+}
+
+template class NSVE<float>;
+template class NSVE<double>;
+

bfps/cpp/full_code/NSVE.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 NSVE_HPP
+#define NSVE_HPP
+
+
+
+#include <cstdlib>
+#include "base.hpp"
+#include "vorticity_equation.hpp"
+#include "full_code/direct_numerical_simulation.hpp"
+
+template <typename rnumber>
+class NSVE: public direct_numerical_simulation
+{
+    public:
+
+        /* parameters that are read in read_parameters */
+        double dt;
+        double famplitude;
+        double fk0;
+        double fk1;
+        int fmode;
+        char forcing_type[512];
+        int histogram_bins;
+        double max_velocity_estimate;
+        double max_vorticity_estimate;
+        double nu;
+
+        /* other stuff */
+        vorticity_equation<rnumber, FFTW> *fs;
+        field<rnumber, FFTW, THREE> *tmp_vec_field;
+        field<rnumber, FFTW, ONE> *tmp_scal_field;
+
+
+        NSVE(
+                const MPI_Comm COMMUNICATOR,
+                const std::string &simulation_name):
+            direct_numerical_simulation(
+                    COMMUNICATOR,
+                    simulation_name){}
+        ~NSVE(){}
+
+        int initialize(void);
+        int main_loop(void);
+        int finalize(void);
+
+        int read_parameters(void);
+        int do_stats(void);
+};
+
+#endif//NSVE_HPP
+

bfps/cpp/full_code/NSVEp.cpp

+#include <string>
+#include <cmath>
+#include "NSVEp.hpp"
+
+
+template <typename rnumber>
+int NSVEp<rnumber>::initialize(void)
+{
+    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);
+        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->fs = new vorticity_equation<rnumber, FFTW>(
+            simname.c_str(),
+            nx, ny, nz,
+            dkx, dky, dkz,
+            FFTW_MEASURE);
+    this->tmp_vec_field = new field<rnumber, FFTW, THREE>(
+            nx, ny, nz,
+            this->comm,
+            FFTW_MEASURE);
+
+
+    this->fs->checkpoints_per_file = checkpoints_per_file;
+    this->fs->nu = nu;
+    this->fs->fmode = fmode;
+    this->fs->famplitude = famplitude;
+    this->fs->fk0 = fk0;
+    this->fs->fk1 = fk1;
+    strncpy(this->fs->forcing_type, forcing_type, 128);
+    this->fs->iteration = this->iteration;
+    this->fs->checkpoint = this->checkpoint;
+
+    this->fs->cvorticity->real_space_representation = false;
+    this->fs->io_checkpoint();
+
+    if (this->myrank == 0 && this->iteration == 0)
+        this->fs->kk->store(stat_file);
+
+    this->ps = particles_system_builder(
+                fs->cvelocity,              // (field object)
+                fs->kk,                     // (kspace object, contains dkx, dky, dkz)
+                tracers0_integration_steps, // to check coherency between parameters and hdf input file (nb rhs)
+                nparticles,                 // to check coherency between parameters and hdf input file
+                fs->get_current_fname(),    // particles input filename
+                std::string("/tracers0/state/") + std::to_string(fs->iteration), // dataset name for initial input
+                std::string("/tracers0/rhs/")  + std::to_string(fs->iteration), // dataset name for initial input
+                tracers0_neighbours,        // parameter (interpolation no neighbours)
+                tracers0_smoothness,        // parameter
+                this->comm,
+                fs->iteration+1);
+    this->particles_output_writer_mpi = new particles_output_hdf5<double,3,3>(
+                MPI_COMM_WORLD,
+                "tracers0",
+                nparticles,
+                tracers0_integration_steps);
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVEp<rnumber>::main_loop(void)
+{
+    clock_t time0, time1;
+    double time_difference, local_time_difference;
+    time0 = clock();
+    bool stop_code_now = false;
+    int max_iter = (this->iteration + this->niter_todo -
+                    (this->iteration % this->niter_todo));
+    for (; this->iteration < max_iter;)
+    {
+    #ifdef USE_TIMINGOUTPUT
+        const std::string loopLabel = ("code::main_start::loop-" +
+                                       std::to_string(this->iteration));
+        TIMEZONE(loopLabel.c_str());
+    #endif
+        if (this->iteration % this->niter_stat == 0)
+            this->do_stats();
+
+        this->fs->compute_velocity(fs->cvorticity);
+        this->fs->cvelocity->ift();
+        this->ps->completeLoop(dt);
+
+        this->fs->step(dt);
+        this->iteration = this->fs->iteration;
+        if (this->fs->iteration % this->niter_out == 0)
+        {
+            this->fs->io_checkpoint(false);
+            this->particles_output_writer_mpi->open_file(fs->get_current_fname());
+            this->particles_output_writer_mpi->save(
+                    this->ps->getParticlesPositions(),
+                    this->ps->getParticlesRhs(),
+                    this->ps->getParticlesIndexes(),
+                    this->ps->getLocalNbParticles(),
+                    this->fs->iteration);
+            this->particles_output_writer_mpi->close_file();
+            this->checkpoint = this->fs->checkpoint;
+            this->write_iteration();
+        }
+        if (stop_code_now)
+            break;
+        time1 = clock();
+        local_time_difference = ((
+                (unsigned int)(time1 - time0)) /
+                ((double)CLOCKS_PER_SEC));
+        time_difference = 0.0;
+        MPI_Allreduce(
+                &local_time_difference,
+                &time_difference,
+                1,
+                MPI_DOUBLE,
+                MPI_SUM,
+                MPI_COMM_WORLD);
+        if (this->myrank == 0)
+            std::cout << "iteration " << iteration << " took " << time_difference/nprocs << " seconds" << std::endl;
+        if (this->myrank == 0)
+            std::cerr << "iteration " << iteration << " took " << time_difference/nprocs << " seconds" << std::endl;
+        time0 = time1;
+    }
+    if (this->iteration % this->niter_stat == 0)
+        this->do_stats();
+    time1 = clock();
+    local_time_difference = ((
+            (unsigned int)(time1 - time0)) /
+            ((double)CLOCKS_PER_SEC));
+    time_difference = 0.0;
+    MPI_Allreduce(
+            &local_time_difference,
+            &time_difference,
+            1,
+            MPI_DOUBLE,
+            MPI_SUM,
+            MPI_COMM_WORLD);
+    if (this->myrank == 0)
+        std::cout << "iteration " << iteration << " took " << time_difference/nprocs << " seconds" << std::endl;
+    if (this->myrank == 0)
+        std::cerr << "iteration " << iteration << " took " << time_difference/nprocs << " seconds" << std::endl;
+    if (this->iteration % this->niter_out != 0)
+    {
+        this->fs->io_checkpoint(false);
+        this->particles_output_writer_mpi->open_file(fs->get_current_fname());
+        this->particles_output_writer_mpi->save(
+                this->ps->getParticlesPositions(),
+                this->ps->getParticlesRhs(),
+                this->ps->getParticlesIndexes(),
+                this->ps->getLocalNbParticles(),
+                this->fs->iteration);
+        this->particles_output_writer_mpi->close_file();
+        this->checkpoint = this->fs->checkpoint;
+        this->write_iteration();
+    }
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVEp<rnumber>::finalize(void)
+{
+    if (this->myrank == 0)
+        H5Fclose(this->stat_file);
+    delete this->fs;
+    delete this->tmp_vec_field;
+    this->ps.release();
+    delete this->particles_output_writer_mpi;
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVEp<rnumber>::do_stats()
+{
+    hid_t stat_group;
+    if (this->myrank == 0)
+        stat_group = H5Gopen(
+                this->stat_file,
+                "statistics",
+                H5P_DEFAULT);
+    else
+        stat_group = 0;
+    fs->compute_velocity(fs->cvorticity);
+    *tmp_vec_field = fs->cvelocity->get_cdata();
+    tmp_vec_field->compute_stats(
+            fs->kk,
+            stat_group,
+            "velocity",
+            fs->iteration / niter_stat,
+            max_velocity_estimate/sqrt(3));
+
+    *tmp_vec_field = fs->cvorticity->get_cdata();
+    tmp_vec_field->compute_stats(
+            fs->kk,
+            stat_group,
+            "vorticity",
+            fs->iteration / niter_stat,
+            max_vorticity_estimate/sqrt(3));
+
+    if (this->myrank == 0)
+        H5Gclose(stat_group);
+
+    if (myrank == 0)
+    {
+        std::string fname = (
+                std::string("stop_") +
+                std::string(simname));
+        {
+            struct stat file_buffer;
+            this->stop_code_now = (
+                    stat(fname.c_str(), &file_buffer) == 0);
+        }
+    }
+    MPI_Bcast(
+            &this->stop_code_now,
+            1,
+            MPI_C_BOOL,
+            0,
+            MPI_COMM_WORLD);
+    return EXIT_SUCCESS;
+}
+
+template class NSVEp<float>;
+template class NSVEp<double>;
+

bfps/cpp/full_code/NSVEp.hpp

+/**********************************************************************
+*                                                                     *
+*  Copyright 2017 Max Planck Institute                                *
+*                 for Dynamics and Self-Organization                  *
+*                                                                     *
+*  This file is part of bfps.                                         *
+*                                                                     *