Merge branch 'develop'

bfps/cpp/interpolator.hpp

@@ -27,23 +27,12 @@
 #include "field_descriptor.hpp"
 #include "fftw_tools.hpp"
 #include "fluid_solver_base.hpp"
-#include "spline_n1.hpp"
-#include "spline_n2.hpp"
-#include "spline_n3.hpp"
-#include "spline_n4.hpp"
-#include "spline_n5.hpp"
-#include "spline_n6.hpp"
-#include "Lagrange_polys.hpp"
+#include "interpolator_base.hpp"
 
 #ifndef INTERPOLATOR
 
 #define INTERPOLATOR
 
-typedef void (*base_polynomial_values)(
-        const int derivative,
-        const double fraction,
-        double *__restrict__ destination);
-
 template <class rnumber, int interp_neighbours>
 class interpolator
 {

bfps/cpp/interpolator_base.hpp

+/**********************************************************************
+*                                                                     *
+*  Copyright 2015 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 "spline_n1.hpp"
+#include "spline_n2.hpp"
+#include "spline_n3.hpp"
+#include "spline_n4.hpp"
+#include "spline_n5.hpp"
+#include "spline_n6.hpp"
+#include "Lagrange_polys.hpp"
+
+#ifndef INTERPOLATOR_BASE
+
+#define INTERPOLATOR_BASE
+
+typedef void (*base_polynomial_values)(
+        const int derivative,
+        const double fraction,
+        double *__restrict__ destination);
+
+#endif//INTERPOLATOR_BASE
+

bfps/cpp/particles.hpp

@@ -29,6 +29,7 @@
 #include <iostream>
 #include <hdf5.h>
 #include "base.hpp"
+#include "particles_base.hpp"
 #include "fluid_solver_base.hpp"
 #include "interpolator.hpp"
 
@@ -36,9 +37,6 @@
 
 #define PARTICLES
 
-/* particle types */
-enum particle_types {VELOCITY_TRACER};
-
 template <int particle_type, class rnumber, bool multistep, int interp_neighbours>
 class particles
 {

bfps/cpp/particles_base.hpp

+/**********************************************************************
+*                                                                     *
+*  Copyright 2015 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 "interpolator_base.hpp"
+
+#ifndef PARTICLES_BASE
+
+#define PARTICLES_BASE
+
+/* particle types */
+enum particle_types {VELOCITY_TRACER};
+
+#endif//PARTICLES_BASE
+

bfps/cpp/rFFTW_interpolator.cpp

@@ -32,15 +32,13 @@
 template <class rnumber, int interp_neighbours>
 rFFTW_interpolator<rnumber, interp_neighbours>::rFFTW_interpolator(
         fluid_solver_base<rnumber> *fs,
-        base_polynomial_values BETA_POLYS)
+        base_polynomial_values BETA_POLYS,
+        rnumber *FIELD)
 {
     this->descriptor = fs->rd;
     this->field_size = 2*fs->cd->local_size;
     this->compute_beta = BETA_POLYS;
-    if (sizeof(rnumber) == 4)
-        this->field = (rnumber*)((void*)fftwf_alloc_real(this->field_size));
-    else if (sizeof(rnumber) == 8)
-        this->field = (rnumber*)((void*)fftw_alloc_real(this->field_size));
+    this->field = FIELD;
 
     // compute dx, dy, dz;
     this->dx = 4*acos(0) / (fs->dkx*this->descriptor->sizes[2]);
@@ -59,24 +57,9 @@ rFFTW_interpolator<rnumber, interp_neighbours>::rFFTW_interpolator(
 template <class rnumber, int interp_neighbours>
 rFFTW_interpolator<rnumber, interp_neighbours>::~rFFTW_interpolator()
 {
-    if (sizeof(rnumber) == 4)
-        fftwf_free((float*)((void*)this->field));
-    else if (sizeof(rnumber) == 8)
-        fftw_free((double*)((void*)this->field));
     delete[] this->compute;
 }
 
-template <class rnumber, int interp_neighbours>
-int rFFTW_interpolator<rnumber, interp_neighbours>::read_rFFTW(void *void_src)
-{
-    rnumber *src = (rnumber*)void_src;
-    /* do big copy of middle stuff */
-    std::copy(src,
-              src + this->field_size,
-              this->field);
-    return EXIT_SUCCESS;
-}
-
 template <class rnumber, int interp_neighbours>
 void rFFTW_interpolator<rnumber, interp_neighbours>::get_grid_coordinates(
         const int nparticles,

bfps/cpp/rFFTW_interpolator.hpp

@@ -27,14 +27,7 @@
 #include "field_descriptor.hpp"
 #include "fftw_tools.hpp"
 #include "fluid_solver_base.hpp"
-#include "spline_n1.hpp"
-#include "spline_n2.hpp"
-#include "spline_n3.hpp"
-#include "spline_n4.hpp"
-#include "spline_n5.hpp"
-#include "spline_n6.hpp"
-#include "Lagrange_polys.hpp"
-#include "interpolator.hpp"
+#include "interpolator_base.hpp"
 
 #ifndef RFFTW_INTERPOLATOR
 
@@ -67,7 +60,8 @@ class rFFTW_interpolator
 
         rFFTW_interpolator(
                 fluid_solver_base<rnumber> *FSOLVER,
-                base_polynomial_values BETA_POLYS);
+                base_polynomial_values BETA_POLYS,
+                rnumber *FIELD_DATA);
         ~rFFTW_interpolator();
 
         /* map real locations to grid coordinates */
@@ -90,7 +84,6 @@ class rFFTW_interpolator
                 const double *__restrict__ xx,
                 double *__restrict__ dest,
                 const int *deriv = NULL);
-        int read_rFFTW(void *src);
 };
 
 #endif//RFFTW_INTERPOLATOR

bfps/cpp/rFFTW_particles.hpp

@@ -29,9 +29,9 @@
 #include <iostream>
 #include <hdf5.h>
 #include "base.hpp"
+#include "particles_base.hpp"
 #include "fluid_solver_base.hpp"
 #include "rFFTW_interpolator.hpp"
-#include "particles.hpp"
 
 #ifndef RFFTW_PARTICLES
 

bfps/fluid_base.py

@@ -130,9 +130,10 @@ class fluid_particle_base(bfps.code):
                              'hsize_t dset;\n')
         for key in ['state', 'velocity', 'acceleration']:
             self.definitions += ('full_name = (std::string(name) + std::string("/{0}"));\n'.format(key) +
+                                 'if (H5Lexists(g_id, full_name.c_str(), H5P_DEFAULT))\n{\n' +
                                  'dset = H5Dopen(g_id, full_name.c_str(), H5P_DEFAULT);\n' +
                                  'grow_single_dataset(dset, niter_todo/niter_part);\n' +
-                                 'H5Dclose(dset);\n')
+                                 'H5Dclose(dset);\n}\n')
         self.definitions += ('full_name = (std::string(name) + std::string("/rhs"));\n' +
                              'if (H5Lexists(g_id, full_name.c_str(), H5P_DEFAULT))\n{\n' +
                              'dset = H5Dopen(g_id, full_name.c_str(), H5P_DEFAULT);\n' +
@@ -445,25 +446,24 @@ class fluid_particle_base(bfps.code):
                                      dtype = np.int64,
                                      compression = 'gzip')
             for s in range(self.particle_species):
-                if self.parameters['tracers{0}_integration_method'.format(s)] == 'AdamsBashforth':
-                    time_chunk = 2**20 // (8*3*
-                                           self.parameters['nparticles']*
-                                           self.parameters['tracers{0}_integration_steps'.format(s)])
-                    time_chunk = max(time_chunk, 1)
-                    ofile.create_dataset('particles/tracers{0}/rhs'.format(s),
-                                         (1,
-                                          self.parameters['tracers{0}_integration_steps'.format(s)],
-                                          self.parameters['nparticles'],
-                                          3),
-                                         maxshape = (None,
-                                                     self.parameters['tracers{0}_integration_steps'.format(s)],
-                                                     self.parameters['nparticles'],
-                                                     3),
-                                         chunks =  (time_chunk,
-                                                    self.parameters['tracers{0}_integration_steps'.format(s)],
-                                                    self.parameters['nparticles'],
-                                                    3),
-                                         dtype = np.float64)
+                time_chunk = 2**20 // (8*3*
+                                       self.parameters['nparticles']*
+                                       self.parameters['tracers{0}_integration_steps'.format(s)])
+                time_chunk = max(time_chunk, 1)
+                ofile.create_dataset('particles/tracers{0}/rhs'.format(s),
+                                     (1,
+                                      self.parameters['tracers{0}_integration_steps'.format(s)],
+                                      self.parameters['nparticles'],
+                                      3),
+                                     maxshape = (None,
+                                                 self.parameters['tracers{0}_integration_steps'.format(s)],
+                                                 self.parameters['nparticles'],
+                                                 3),
+                                     chunks =  (time_chunk,
+                                                self.parameters['tracers{0}_integration_steps'.format(s)],
+                                                self.parameters['nparticles'],
+                                                3),
+                                     dtype = np.float64)
                 time_chunk = 2**20 // (8*3*self.parameters['nparticles'])
                 time_chunk = max(time_chunk, 1)
                 ofile.create_dataset(
@@ -474,14 +474,15 @@ class fluid_particle_base(bfps.code):
                     chunks = (time_chunk, self.parameters['nparticles'], 3),
                     maxshape = (None, self.parameters['nparticles'], 3),
                     dtype = np.float64)
-                ofile.create_dataset(
-                    '/particles/tracers{0}/acceleration'.format(s),
-                    (1,
-                     self.parameters['nparticles'],
-                     3),
-                    chunks = (time_chunk, self.parameters['nparticles'], 3),
-                    maxshape = (None, self.parameters['nparticles'], 3),
-                    dtype = np.float64)
+                if self.parameters['tracers{0}_acc_on'.format(s)]:
+                    ofile.create_dataset(
+                        '/particles/tracers{0}/acceleration'.format(s),
+                        (1,
+                         self.parameters['nparticles'],
+                         3),
+                        chunks = (time_chunk, self.parameters['nparticles'], 3),
+                        maxshape = (None, self.parameters['nparticles'], 3),
+                        dtype = np.float64)
             ofile.close()
         return None
 

done.txt

@@ -3,3 +3,10 @@ x 2015-12-23 decide on versioning system
 x 2015-12-24 move get_grid coords to interpolator                    @optimization +v1.0
 x 2015-12-25 get rid of temporal interpolation                       @optimization +v1.0
 x 2015-12-26 call interpolation only when needed                     @optimization +v1.0
+x 2015-12-26 clean up tox files, make sure all tests run             @tests        +v1.0
+x 2016-01-03 check divfree function
+x 2016-01-03 compute kMeta(t) as well
+x 2016-01-03 split library into core and extra                                       @optimization +v1.0
+x 2016-01-07 FFTW interpolator doesn't need its own field                            @optimization +v1.0 +particle_api
+x 2016-01-08 simplify tracer/field addition mechanism                                @design +v1.0 +particle_api
+x 2016-01-08 add stat choice parameter to add_particles                              @design +v1.0 +particle_api

setup.py

@@ -68,9 +68,7 @@ print(VERSION)
 src_file_list = ['field_descriptor',
                  'fluid_solver_base',
                  'fluid_solver',
-                 'interpolator',
                  'rFFTW_interpolator',
-                 'particles',
                  'rFFTW_particles',
                  'fftw_tools',
                  'spline_n1',
@@ -81,7 +79,11 @@ src_file_list = ['field_descriptor',
                  'spline_n6',
                  'Lagrange_polys']
 
-header_list = ['cpp/base.hpp'] + ['cpp/' + fname + '.hpp' for fname in src_file_list]
+header_list = (['cpp/base.hpp',
+                'cpp/particles_base.hpp',
+                'cpp/interpolator_base.hpp'] +
+               ['cpp/' + fname + '.hpp'
+                for fname in src_file_list])
 
 with open('MANIFEST.in', 'w') as manifest_in_file:
     for fname in ['bfps/cpp/' + fname + '.cpp' for fname in src_file_list] + header_list:

tests/base.py

@@ -102,29 +102,19 @@ def launch(
     c.parameters['famplitude'] = 0.2
     c.fill_up_fluid_code()
     if c.parameters['nparticles'] > 0:
-        c.add_particle_fields(
-                name = 'regular',
+        c.add_3D_rFFTW_field(name = 'rFFTW_acc')
+        c.add_interpolator(
+                name = 'spline',
                 neighbours = opt.neighbours,
                 smoothness = opt.smoothness)
-        c.add_particle_fields(kcut = 'fs->kM/2', name = 'filtered', neighbours = opt.neighbours)
         c.add_particles(
-                kcut = 'fs->kM/2',
-                integration_steps = 1,
-                fields_name = 'filtered')
-        #for integr_steps in range(1, 7):
-        #    c.add_particles(
-        #            integration_steps = integr_steps,
-        #            neighbours = opt.neighbours,
-        #            smoothness = opt.smoothness,
-        #            fields_name = 'regular')
-        for info in [(2, 'AdamsBashforth'),
-                     (3, 'AdamsBashforth'),
-                     (4, 'AdamsBashforth'),
-                     (6, 'AdamsBashforth')]:
-            c.add_particles(
-                    integration_steps = info[0],
-                    integration_method = info[1],
-                    fields_name = 'regular')
+                kcut = ['fs->kM/2', 'fs->kM/3'],
+                integration_steps = 3,
+                interpolator = 'spline')
+        c.add_particles(
+                integration_steps = [2, 3, 4, 6],
+                interpolator = 'spline',
+                acc_name = 'rFFTW_acc')
     c.finalize_code()
     c.write_src()
     c.write_par()
@@ -158,6 +148,8 @@ def launch(
 
 
 def acceleration_test(c, m = 3, species = 0):
+    if not c.parameters['tracers{0}_acc_on'.format(species)]:
+        return None
     import numpy as np
     import matplotlib.pyplot as plt
     from bfps.tools import get_fornberg_coeffs
@@ -182,13 +174,12 @@ def acceleration_test(c, m = 3, species = 0):
     pid = np.argmin(SNR(num_acc1, acc[n+1:-n-1]))
     pars = d['parameters']
     to_print = (
-            'integration={0}, steps={1}, interp={2}, neighbours={3}, '.format(
-                pars['tracers{0}_integration_method'.format(species)].value,
+            'steps={0}, interp={1}, neighbours={2}, '.format(
                 pars['tracers{0}_integration_steps'.format(species)].value,
-                pars[str(pars['tracers{0}_field'.format(species)].value) + '_type'].value,
-                pars[str(pars['tracers{0}_field'.format(species)].value) + '_neighbours'].value))
-    if 'spline' in pars['tracers{0}_field'.format(species)].value:
-        to_print += 'smoothness = {0}, '.format(pars[str(pars['tracers{0}_field'.format(species)].value) + '_smoothness'].value)
+                pars[str(pars['tracers{0}_interpolator'.format(species)].value) + '_type'].value,
+                pars[str(pars['tracers{0}_interpolator'.format(species)].value) + '_neighbours'].value))
+    if 'spline' in pars['tracers{0}_interpolator'.format(species)].value:
+        to_print += 'smoothness = {0}, '.format(pars[str(pars['tracers{0}_interpolator'.format(species)].value) + '_smoothness'].value)
     to_print += (
             'SNR d1p-vel={0:.3f}, d1v-acc={1:.3f}, d2p-acc={2:.3f}'.format(
                 np.mean(SNR(num_vel1, vel[n+1:-n-1])),

tests/test_convergence.py

-#! /usr/bin/env python2
 #######################################################################
 #                                                                     #
 #  Copyright 2015 Max Planck Institute                                #
@@ -30,7 +29,7 @@ import h5py
 from mpl_toolkits.mplot3d import axes3d
 import matplotlib.pyplot as plt
 
-from test_base import *
+from base import *
 
 def convergence_test(
         opt,
@@ -224,6 +223,15 @@ def convergence_test(
     return c0, c1, c2
 
 if __name__ == '__main__':
-    opt = parser.parse_args()
+    opt = parser.parse_args(
+            ['-n', '16',
+             '--run',
+             '--initialize',
+             '--ncpu', '2',
+             '--nparticles', '1000',
+             '--niter_todo', '32',
+             '--precision', 'single',
+             '--wd', 'data/single'] +
+            sys.argv[1:])
     convergence_test(opt, launch)
 

tests/test_frozen_field.py

-#! /usr/bin/env python2
 #######################################################################
 #                                                                     #
 #  Copyright 2015 Max Planck Institute                                #
@@ -25,7 +24,7 @@
 
 
 
-from test_base import *
+from base import *
 
 class FrozenFieldParticles(bfps.NavierStokes):
     def __init__(
@@ -33,12 +32,16 @@ class FrozenFieldParticles(bfps.NavierStokes):
             name = 'FrozenFieldParticles',
             work_dir = './',
             simname = 'test',
-            fluid_precision = 'single'):
+            frozen_fields = True,
+            fluid_precision = 'single',
+            use_fftw_wisdom = False):
         super(FrozenFieldParticles, self).__init__(
                 name = name,
                 work_dir = work_dir,
                 simname = simname,
-                fluid_precision = fluid_precision)
+                fluid_precision = fluid_precision,
+                frozen_fields = frozen_fields,
+                use_fftw_wisdom = use_fftw_wisdom)
     def fill_up_fluid_code(self):
         self.fluid_includes += '#include <cstring>\n'
         self.fluid_variables += 'fluid_solver<{0}> *fs;\n'.format(self.C_dtype)

tests/test_interpolation.py

-#! /usr/bin/env python2
 #######################################################################
 #                                                                     #
 #  Copyright 2015 Max Planck Institute                                #
@@ -26,35 +25,44 @@
 
 
 import numpy as np
-from test_base import *
+from base import *
 import matplotlib.pyplot as plt
 
-class test_interpolation(bfps.NavierStokes):
-    def __init__(
-            self,
-            name = 'test_interpolation',
-            work_dir = './',
-            simname = 'test'):
-        super(test_interpolation, self).__init__(
-                work_dir = work_dir,
-                simname = simname,
-                name = name,
-                frozen_fields = True)
-        return None
-
 if __name__ == '__main__':
-    opt = parser.parse_args()
-    c = test_interpolation(work_dir = opt.work_dir + '/io')
+    opt = parser.parse_args(
+            ['-n', '32',
+             '--run',
+             '--ncpu', '2',
+             '--initialize',
+             '--neighbours', '4',
+             '--smoothness', '1',
+             '--nparticles', '4225',
+             '--niter_todo', '1',
+             '--niter_stat', '1',
+             '--niter_part', '1',
+             '--niter_out', '1',
+             '--wd', 'interp'] +
+            sys.argv[1:])
+    c = bfps.NavierStokes(
+            name = 'test_interpolation',
+            use_fftw_wisdom = False,
+            frozen_fields = True)
     c.pars_from_namespace(opt)
+    c.fill_up_fluid_code()
+    c.add_particle_fields(
+            name = 'n{0}m{1}'.format(opt.neighbours, opt.smoothness),
+            neighbours = opt.neighbours,
+            smoothness = opt.smoothness)
+    c.add_particle_fields(
+            name = 'n{0}'.format(opt.neighbours),
+            neighbours = opt.neighbours,
+            interp_type = 'Lagrange')
     c.add_particles(
             integration_steps = 1,
-            neighbours = 4,
-            smoothness = 1)
+            fields_name = 'n{0}m{1}'.format(opt.neighbours, opt.smoothness))
     c.add_particles(
             integration_steps = 1,
-            neighbours = 4,
-            interp_type = 'Lagrange')
-    c.fill_up_fluid_code()
+            fields_name = 'n{0}'.format(opt.neighbours))
     c.finalize_code()
     c.write_src()
     c.write_par()

tests/test_io.py

-#! /usr/bin/env python2
 #######################################################################
 #                                                                     #
 #  Copyright 2015 Max Planck Institute                                #
@@ -25,7 +24,7 @@
 
 
 
-from test_base import *
+from base import *
 
 class test_io(bfps.code):
     def __init__(
@@ -43,7 +42,12 @@ class test_io(bfps.code):
         return None
 
 if __name__ == '__main__':
-    opt = parser.parse_args()
+    opt = parser.parse_args(
+            ['-n', '32',
+             '--run',
+             '--initialize',
+             '--ncpu', '2'] +
+            sys.argv[1:])
     c = test_io(work_dir = opt.work_dir + '/io')
     c.write_src()
     c.write_par()

tests/test_particles.py

@@ -28,7 +28,7 @@
 import numpy as np
 import matplotlib.pyplot as plt
 
-from test_base import *
+from base import *
 
 from test_frozen_field import FrozenFieldParticles
 from test_convergence import convergence_test
@@ -147,7 +147,17 @@ class err_finder:
         return errlist
 
 if __name__ == '__main__':
-    opt = parser.parse_args()
+    opt = parser.parse_args(
+            ['-n', '16',
+             '--run',
+             '--initialize',
+             '--frozen',
+             '--ncpu', '2',
+             '--nparticles', '128',
+             '--niter_todo', '16',
+             '--precision', 'single',
+             '--wd', 'data/single'] +
+            sys.argv[1:])
     if opt.precision == 'single':
         dtype = np.complex64
     elif opt.precision == 'double':

tests/test_plain.py

@@ -65,6 +65,7 @@ if __name__ == '__main__':
     opt = parser.parse_args(
             ['-n', '32',
              '--run',
+             '--initialize',
              '--ncpu', '2',
              '--nparticles', '1000',
              '--niter_todo', '48',

tests/test_time_step.py

-#! /usr/bin/env python2
 #######################################################################
 #                                                                     #
 #  Copyright 2015 Max Planck Institute                                #
@@ -30,7 +29,7 @@ import h5py
 from mpl_toolkits.mplot3d import axes3d
 import matplotlib.pyplot as plt
 
-from test_base import *
+from base import *
 
 def convergence_test(
         opt,
@@ -93,10 +92,19 @@ def convergence_test(