Merge branch 'feature/projects' into develop

projects/sedimenting_rods/GCDNS.py

+import os
+import sys
+import argparse
+import getpass
+from packaging import version
+
+username = getpass.getuser()
+homefolder = os.path.expanduser('~')
+
+here = os.path.normcase(__file__)
+
+cpp_location = os.path.dirname(here)
+
+import numpy as np
+import h5py
+
+import TurTLE
+#assert(version.parse(TurTLE.__version__) > version.parse('3.16'))
+import TurTLE._base
+import TurTLE.DNS
+from TurTLE._code import _code
+
+
+class GCDNS(TurTLE.DNS):
+    def __init__(
+            self,
+            **kwargs):
+        self.dns_type = 'many_NSVEcylinder_collisions'
+        TurTLE.DNS.__init__(self, **kwargs)
+        self.extra_parameters = {} #Added this here to get it to work
+        self.extra_parameters[self.dns_type] = {
+            'niter_part' : int(10),
+            'niter_part_fine_period' : int(10),
+            'niter_part_fine_duration' : int(0),
+            'nparticles' : np.array([100],dtype=int),
+            'tracers0_integration_steps' : int(2),
+            'tracers0_neighbours' : int(3),
+            'tracers0_smoothness' : int(2),
+            'tracers0_enable_p2p' : int(1),
+            'tracers0_enable_inner' : int(0),
+            'tracers0_enable_vorticity_omega' : int(0),
+            'tracers0_cutoff' : np.array([0.5],dtype='float64'),
+            'tracers0_lambda' : np.array([100.0],dtype='float64'),
+            'tracers0_beta0' : np.array([0.0],dtype='float64'),
+            'tracers0_beta1' : np.array([0.0],dtype='float64'),
+            'tracers0_D_rho' : np.array([1.0],dtype='float64'),
+            'tracers0_u_sed' : np.array([1.0],dtype='float64'),
+            'l_factor' : float(1.0),
+            't_factor' : float(1.0),
+            'many_integration_steps' : np.array([2]).astype(np.int),
+                }
+        return None
+    def write_src(
+            self):
+        self.version_message = (
+                '/***********************************************************************\n' +
+                '* this code automatically generated by TurTLE\n' +
+                '* version {0}\n'.format(TurTLE.__version__) +
+                '***********************************************************************/\n\n\n')
+        self.include_list = [
+                '"base.hpp"',
+                '"scope_timer.hpp"',
+                '"fftw_interface.hpp"',
+                '"full_code/main_code.hpp"',
+                '"full_code/NSVEparticles.hpp"',
+                '<cmath>',
+                '<iostream>',
+                '<hdf5.h>',
+                '<string>',
+                '<cstring>',
+                '<fftw3-mpi.h>',
+                '<omp.h>',
+                '<cfenv>',
+                '<cstdlib>']
+        self.main = """
+            int main(int argc, char *argv[])
+            {{
+                bool fpe = (
+                    (getenv("BFPS_FPE_OFF") == nullptr) ||
+                    (getenv("BFPS_FPE_OFF") != std::string("TRUE")));
+                return main_code< {0} >(argc, argv, fpe);
+            }}
+            """.format(self.dns_type + '<{0}>'.format(self.C_field_dtype))
+        self.includes = '\n'.join(
+                ['#include ' + hh
+                 for hh in self.include_list])
+        self.definitions += open(
+                os.path.join(
+                    cpp_location, 'particles_inner_computer_sedimenting_rods.hpp'), 'r').read()
+        self.definitions += open(
+                os.path.join(
+                    cpp_location, 'particles_inner_computer_sedimenting_rods.cpp'), 'r').read()
+        self.definitions += open(
+            os.path.join(
+                cpp_location, 'many_NSVEcylinder_collisions.hpp'), 'r').read()
+        self.definitions += open(
+            os.path.join(
+                cpp_location, 'many_NSVEcylinder_collisions.cpp'), 'r').read()
+        with open(self.name + '.cpp', 'w') as outfile:
+            outfile.write(self.version_message + '\n\n')
+            outfile.write(self.includes + '\n\n')
+            outfile.write(self.definitions + '\n\n')
+            outfile.write(self.main + '\n')
+        self.check_current_vorticity_exists = True
+        return None
+
+    def generate_tracer_state(
+            self,
+            rseed = None,
+            species = 0,
+            integration_steps = None,
+            ncomponents = 3
+            ):
+        try:
+            if type(integration_steps) == type(None):
+                integration_steps = self.NSVEp_extra_parameters['tracers0_integration_steps']
+            if 'tracers{0}_integration_steps'.format(species) in self.parameters.keys():
+                integration_steps = self.parameters['tracers{0}_integration_steps'.format(species)]
+            if self.dns_type == 'static_field_with_cylinder_collisions' or self.dns_type == 'NSVEcylinder_collisions' or self.dns_type == 'many_NSVEcylinder_collisions':
+                ncomponents = 6
+            with h5py.File(self.get_checkpoint_0_fname(), 'a') as data_file:
+                nn = self.parameters['nparticles'][species]
+                if not 'tracers{0}'.format(species) in data_file.keys():
+                    data_file.create_group('tracers{0}'.format(species))
+                    data_file.create_group('tracers{0}/rhs'.format(species))
+                    data_file.create_group('tracers{0}/state'.format(species))
+                data_file['tracers{0}/rhs'.format(species)].create_dataset(
+                        '0',
+                        shape = (integration_steps, nn, ncomponents,),
+                        dtype = np.float)
+                dset = data_file['tracers{0}/state'.format(species)].create_dataset(
+                        '0',
+                        shape = (nn, ncomponents,),
+                        dtype = np.float)
+                if not type(rseed) == type(None):
+                    np.random.seed(rseed)
+                cc = int(0)
+                batch_size = int(1e6)
+                def get_random_phases(npoints):
+                    return np.random.random(
+                                    (npoints, 3))*2*np.pi
+                def get_random_versors(npoints):
+                    if self.parallel_rods:
+                        bla = np.zeros((npoints, 3),dtype='float64')
+                        bla[:,1]=1.0
+                    else:
+                        bla = np.random.normal(
+                                size = (npoints, 3))
+                        bla  /= np.sum(bla**2, axis = 1)[:, None]**.5
+                    return bla
+                while nn > 0:
+                    if nn > batch_size:
+                        dset[cc*batch_size:(cc+1)*batch_size, :3] = get_random_phases(batch_size)
+                        if dset.shape[1] == 6:
+                            dset[cc*batch_size:(cc+1)*batch_size, 3:] = get_random_versors(batch_size)
+                        nn -= batch_size
+                    else:
+                        dset[cc*batch_size:cc*batch_size+nn, :3] = get_random_phases(nn)
+                        if dset.shape[1] == 6:
+                            dset[cc*batch_size:cc*batch_size+nn, 3:] = get_random_versors(nn)
+                        nn = 0
+                    cc += 1
+        except Exception as e:
+            print('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!WARNING')
+            print('WARNING: when generating particles, caught the exception ', e)
+            print('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!WARNING')
+        return None
+
+    def launch(
+            self,
+            args = [],
+            test_collisions = False,
+            parallel_rods = False,
+            **kwargs
+            ):
+        self.parallel_rods = parallel_rods
+        opt = self.prepare_launch(args = args)
+        if not os.path.exists(self.get_data_file_name()):
+            self.write_par()
+            self.generate_initial_condition(opt = opt)
+            for species in range(self.parameters['many_integration_steps'].shape[0]):
+                self.generate_tracer_state(
+                        species = species,
+                        ncomponents = 6,
+                        integration_steps = self.parameters['many_integration_steps'][species-1]
+                        )
+                with h5py.File(self.get_particle_file_name(), 'a') as particle_file:
+                    particle_file.create_group('tracers{0}/position'.format(species))
+                    particle_file.create_group('tracers{0}/orientation'.format(species))
+                    particle_file.create_group('tracers{0}/velocity'.format(species))
+                    particle_file.create_group('tracers{0}/velocity_gradient'.format(species))
+                    particle_file.create_group('tracers{0}/acceleration'.format(species))
+                with h5py.File(os.path.join(self.work_dir, self.simname + '_colliding_particle_subset.h5'), 'a') as subset_file:
+                    subset_file.create_group('subset_tracers{0}/labels'.format(species))
+                    subset_file.create_group('subset_tracers{0}/position'.format(species))
+                    subset_file.create_group('subset_tracers{0}/orientation'.format(species))
+                    subset_file.create_group('subset_tracers{0}/velocity'.format(species))
+                    subset_file.create_group('subset_tracers{0}/velocity_gradient'.format(species))
+                    subset_file.create_group('subset_tracers{0}/acceleration'.format(species))
+        # Create collisions file
+        if not os.path.exists(os.path.join(self.work_dir, self.simname + '_collisions.h5')):
+            with h5py.File(os.path.join(self.work_dir, self.simname + '_collisions.h5'), 'w') as ofile:
+                for i in range(self.parameters['many_integration_steps'].size):
+                    ofile.create_group('/collisions{0}'.format(i))
+                    ofile.create_group('/new_collisions{0}'.format(i))
+                    ofile.create_group('/same_collisions{0}'.format(i))
+                    ofile.create_group('/old_collisions{0}'.format(i))
+        if not os.path.exists(os.path.join(self.work_dir, self.simname + '_collisions_PP.h5')):
+            with h5py.File(os.path.join(self.work_dir, self.simname + '_collisions_PP.h5'), 'w') as ofile:
+                for i in range(self.parameters['many_integration_steps'].size):
+                    ofile.create_group('/new_collisions{0}'.format(i))
+        self.run(
+                nb_processes = opt.nb_processes,
+                nb_threads_per_process = opt.nb_threads_per_process,
+                njobs = opt.njobs,
+                hours = opt.minutes // 60,
+                minutes = opt.minutes % 60,
+                no_submit = opt.no_submit,
+                no_debug = opt.no_debug)
+        return None
+
+def main():
+    bla = GCDNS()
+    bla.launch(['many_NSVEcylinder_collisions'] + sys.argv[1:])
+    return None
+
+if __name__ == '__main__':
+    main()

projects/sedimenting_rods/README.rst

+Collision counter for sedimenting rods
+--------------------------------------
+
+Please see the PNAS publication
+José-Agustı́n Arguedas-Leiva, Jonasz Słomka, Cristian C. Lalescu, Roman Stocker, and Michael Wilczek
+*Elongation enhances encounter rates betweenin turbulence*
+
+In this folder we include the code used to generate the results
+presented in the paper above.

projects/sedimenting_rods/many_NSVEcylinder_collisions.hpp

+/**********************************************************************
+*                                                                     *
+*  Copyright 2020 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 many_NSVEcylinder_collisions_HPP
+#define many_NSVEcylinder_collisions_HPP
+
+#include <cstdlib>
+#include "base.hpp"
+#include "vorticity_equation.hpp"
+#include "full_code/NSVE.hpp"
+#include "particles/particles_system_builder.hpp"
+#include "particles/particles_output_hdf5.hpp"
+#include "particles/particles_sampling.hpp"
+#include "particles/p2p/p2p_computer.hpp"
+#include "particles/p2p/p2p_ghost_collisions.hpp"
+
+/** \brief Navier-Stokes solver that includes cylinders sedimenting in gravity with a ghost collision counter.
+ *
+ *  Child of Navier Stokes vorticity equation solver, this class calls all the
+ *  methods from `NSVE`, and in addition integrates `sedimenting rods
+ *  in the resulting velocity field.
+ *  Sedimenting rods are very small particles with a density offset to the
+ *  fluid flow. They have an orientation and are idealized as very thin ellongated
+ *  ellipsoids, whose orientation evolves accordint to Jeffery's equations.
+ */
+
+template <typename rnumber>
+class many_NSVEcylinder_collisions: public NSVE<rnumber>
+{
+    public:
+
+        /* parameters that are read in read_parameters */
+        std::vector<int> many_integration_steps;
+        int niter_part;
+        int niter_part_fine_period;//NEW
+        int niter_part_fine_duration;//NEW
+        std::vector<int> nparticles;
+        int tracers0_integration_steps;
+        int tracers0_neighbours;
+        int tracers0_smoothness;
+        std::vector<double> tracers0_cutoff;
+        std::vector<double> tracers0_lambda;
+        std::vector<double> tracers0_D_rho;
+        std::vector<double> tracers0_u_sed;
+        std::vector<double> tracers0_beta0;
+        std::vector<double> tracers0_beta1;
+
+        bool enable_p2p;
+        bool enable_inner;
+        bool enable_vorticity_omega;
+
+        /* other stuff */
+        std::vector<std::unique_ptr<abstract_particles_system_with_p2p<long long int, double, p2p_ghost_collisions<double, long long int>>>> many_ps;
+
+        std::vector<particles_output_hdf5<long long int, double,6>*> particles_output_writer_mpi;
+        std::vector<particles_output_sampling_hdf5<long long int, double, double, 3>*> particles_sample_writer_mpi;
+        // field for sampling velocity gradient
+        field<rnumber, FFTW, THREExTHREE> *nabla_u;
+
+        hid_t cols_file;
+        hid_t cols_file_PP;
+        hid_t subset_file;
+
+        many_NSVEcylinder_collisions(
+                const MPI_Comm COMMUNICATOR,
+                const std::string &simulation_name):
+            NSVE<rnumber>(
+                    COMMUNICATOR,
+                    simulation_name),
+                    tracers0_cutoff(0.314),
+                    tracers0_lambda(100.0),
+                    tracers0_D_rho(1.0),
+                    tracers0_u_sed(1.0),
+                    enable_p2p(true),//interaction has to be turned off by hand in p2p_ghost_collision_base
+                    enable_inner(true),
+                    enable_vorticity_omega(true){}//PUT TRUE HERE AS IN NSVE
+        ~many_NSVEcylinder_collisions(){}
+
+        int initialize(void);
+        int step(void);
+        int finalize(void);
+
+        virtual int read_parameters(void);
+        int write_checkpoint(void);
+        int do_stats(void);
+        
+        std::vector <std::vector <long long int>> collision_pairs_old;
+        std::vector <std::vector <long long int>> collision_pairs_new;
+        std::vector <std::vector <long long int>> collision_pairs_same;
+        std::vector <std::vector <long long int>> collision_pairs_time_0;
+        std::vector <std::vector <long long int>> collision_pairs_time_1;
+  
+void reset_collision_lists(int index){
+        this->collision_pairs_old[index].resize(0);
+        this->collision_pairs_new[index].resize(0);
+        this->collision_pairs_same[index].resize(0);
+        }
+        
+};
+
+#endif//many_NSVEcylinder_collisions_HPP

projects/sedimenting_rods/particles_inner_computer_sedimenting_rods.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 "base.hpp"
+#include "particles/particles_utils.hpp"
+#include <cmath>
+
+template <class real_number, class partsize_t>
+template <int size_particle_positions, int size_particle_rhs>
+void particles_inner_computer_sedimenting_rods<real_number, partsize_t>::compute_interaction(
+        const partsize_t nb_particles,
+        const real_number pos_part[],
+        real_number rhs_part[]) const{
+    static_assert(size_particle_positions == 6, "This kernel works only with 6 values for one particle's position");
+    static_assert(size_particle_rhs == 6, "This kernel works only with 6 values per particle's rhs");
+
+    #pragma omp parallel for
+    for(partsize_t idx_part = 0 ; idx_part < nb_particles ; ++idx_part)
+    {
+        //if (fabs(rhs_part[idx_part*size_particle_rhs + IDXC_X]) > 0 ||
+        //    fabs(rhs_part[idx_part*size_particle_rhs + IDXC_Y]) > 0 ||
+        //    fabs(rhs_part[idx_part*size_particle_rhs + IDXC_Z]) > 0){
+        //    DEBUG_MSG("just entered loop, nonzero sampled velocities for particle %ld are %g %g %g\n",
+        //            idx_part,
+        //            rhs_part[idx_part*size_particle_rhs + IDXC_X],
+        //            rhs_part[idx_part*size_particle_rhs + IDXC_Y],
+        //            rhs_part[idx_part*size_particle_rhs + IDXC_Z]);
+        //            }
+		const double px = pos_part[idx_part*size_particle_positions + 3+IDXC_X];
+		const double py = pos_part[idx_part*size_particle_positions + 3+IDXC_Y];
+		const double pz = pos_part[idx_part*size_particle_positions + 3+IDXC_Z];
+
+        rhs_part[idx_part*size_particle_rhs + IDXC_X] -= D_rho * u_sed * beta1 * px * py;
+        rhs_part[idx_part*size_particle_rhs + IDXC_Y] -= D_rho * u_sed * ( beta0 + beta1 * py*py );
+        rhs_part[idx_part*size_particle_rhs + IDXC_Z] -= D_rho * u_sed * beta1 * pz * py;
+        //if (fabs(rhs_part[idx_part*size_particle_rhs + IDXC_X]) > 0 ||
+        //    fabs(rhs_part[idx_part*size_particle_rhs + IDXC_Y]) > 0 ||
+        //    fabs(rhs_part[idx_part*size_particle_rhs + IDXC_Z]) > 0){
+        //    DEBUG_MSG("computed settling velocity for particle %ld are %g %g %g\n",
+        //            idx_part,
+        //            rhs_part[idx_part*size_particle_rhs + IDXC_X],
+        //            rhs_part[idx_part*size_particle_rhs + IDXC_Y],
+        //            rhs_part[idx_part*size_particle_rhs + IDXC_Z]);
+        //            }
+        //if (idx_part == 0)
+        //{
+        //    DEBUG_MSG("particle %ld, moving with %g %g %g, D_rho = %g, u_sed = %g, beta0 = %g, beta1 = %g\n",
+        //            idx_part,
+        //              rhs_part[idx_part*size_particle_rhs + IDXC_X],
+        //              rhs_part[idx_part*size_particle_rhs + IDXC_Y],
+        //              rhs_part[idx_part*size_particle_rhs + IDXC_Z],
+        //              D_rho, u_sed, beta1, beta0);
+        //    DEBUG_MSG("particle %ld, orientation is %g %g %g\n",
+        //            idx_part,
+        //              px,
+        //              py,
+        //              pz);
+        //    DEBUG_MSG("particle %ld, location is %g %g %g\n",
+        //              idx_part,
+        //              pos_part[idx_part*size_particle_positions + IDXC_X],
+        //              pos_part[idx_part*size_particle_positions + IDXC_Y],
+        //              pos_part[idx_part*size_particle_positions + IDXC_Z]);
+        //}
+    }
+}
+    // for given orientation and right-hand-side, recompute right-hand-side such
+    // that it is perpendicular to the current orientation.
+    // this is the job of the Lagrange multiplier terms, hence the
+    // "add_Lagrange_multipliers" name of the method.
+template <>
+template <>
+void particles_inner_computer_sedimenting_rods<double, long long>::add_Lagrange_multipliers<6,6>(
+        const long long nb_particles,
+        const double pos_part[],
+        double rhs_part[]) const{
+
+    //DEBUG_MSG("hello from add_lagrange_multipliers\n");
+
+        #pragma omp parallel for
+        for(long long idx_part = 0 ; idx_part < nb_particles ; ++idx_part){
+            const long long idx0 = idx_part*6 + 3;
+            const long long idx1 = idx_part*6 + 3;
+            // check that orientation is unit vector:
+            double orientation_size = sqrt(
+                    pos_part[idx0+IDXC_X]*pos_part[idx0+IDXC_X] +
+                    pos_part[idx0+IDXC_Y]*pos_part[idx0+IDXC_Y] +
+                    pos_part[idx0+IDXC_Z]*pos_part[idx0+IDXC_Z]);
+            variable_used_only_in_assert(orientation_size);
+            assert(orientation_size > 0.99);
+            assert(orientation_size < 1.01);
+            // I call "rotation" to be the right hand side of the orientation part of the ODE
+            // project rotation on orientation:
+            double projection = (
+                    pos_part[idx0+IDXC_X]*rhs_part[idx1+IDXC_X] +
+                    pos_part[idx0+IDXC_Y]*rhs_part[idx1+IDXC_Y] +
+                    pos_part[idx0+IDXC_Z]*rhs_part[idx1+IDXC_Z]);
+
+            // now remove parallel bit.
+            rhs_part[idx1+IDXC_X] -= pos_part[idx0+IDXC_X]*projection;
+            rhs_part[idx1+IDXC_Y] -= pos_part[idx0+IDXC_Y]*projection;
+            rhs_part[idx1+IDXC_Z] -= pos_part[idx0+IDXC_Z]*projection;
+            //if (idx_part == 0){
+             //   DEBUG_MSG("add_Lagrange_multipliers particles 0 projection is %g\n", projection);
+            //}
+
+            // DEBUG
+            // sanity check, for debugging purposes
+            // compute dot product between orientation and orientation change
+            //double dotproduct = (
+            //        rhs_part[idx1 + IDXC_X]*pos_part[idx0 + IDXC_X] +
+            //        rhs_part[idx1 + IDXC_Y]*pos_part[idx0 + IDXC_Y] +
+            //        rhs_part[idx1 + IDXC_Z]*pos_part[idx0 + IDXC_Z]);
+            //if (dotproduct > 0.1)
+            //{
+            //    DEBUG_MSG("dotproduct = %g, projection = %g\n"
+            //              "pos_part[%d] = %g, pos_part[%d] = %g, pos_part[%d] = %g\n"
+            //              "rhs_part[%d] = %g, rhs_part[%d] = %g, rhs_part[%d] = %g\n",
+            //            dotproduct,
+            //            projection,
+            //            IDXC_X, pos_part[idx0 + IDXC_X],
+            //            IDXC_Y, pos_part[idx0 + IDXC_Y],
+            //            IDXC_Z, pos_part[idx0 + IDXC_Z],
+            //            IDXC_X, rhs_part[idx1 + IDXC_X],
+            //            IDXC_Y, rhs_part[idx1 + IDXC_Y],
+            //            IDXC_Z, rhs_part[idx1 + IDXC_Z]);
+            //    assert(false);
+            //}
+            //assert(dotproduct <= 0.1);
+        }
+    }
+
+template <>
+template <>
+void particles_inner_computer_sedimenting_rods<double, long long>::compute_interaction_with_extra<6,6,3>(
+        const long long nb_particles,
+        const double pos_part[],
+        double rhs_part[],
+        const double rhs_part_extra[]) const{
+    // call plain compute_interaction first
+    //DEBUG_MSG("hello from compute_interaction_with_extra 3\n");
+    compute_interaction<6, 6>(nb_particles, pos_part, rhs_part);
+
+    // now add vorticity term
+    #pragma omp parallel for
+    for(long long idx_part = 0 ; idx_part < nb_particles ; ++idx_part){
+        // Cross product vorticity/orientation
+        rhs_part[idx_part*6 + 3+IDXC_X] += 0.5*(rhs_part_extra[idx_part*3 + IDXC_Y]*pos_part[idx_part*6 + 3+IDXC_Z] -
+                                                rhs_part_extra[idx_part*3 + IDXC_Z]*pos_part[idx_part*6 + 3+IDXC_Y]);
+        rhs_part[idx_part*6 + 3+IDXC_Y] += 0.5*(rhs_part_extra[idx_part*3 + IDXC_Z]*pos_part[idx_part*6 + 3+IDXC_X] -
+                                                rhs_part_extra[idx_part*3 + IDXC_X]*pos_part[idx_part*6 + 3+IDXC_Z]);
+        rhs_part[idx_part*6 + 3+IDXC_Z] += 0.5*(rhs_part_extra[idx_part*3 + IDXC_X]*pos_part[idx_part*6 + 3+IDXC_Y] -
+                                                rhs_part_extra[idx_part*3 + IDXC_Y]*pos_part[idx_part*6 + 3+IDXC_X]);
+    }
+}
+
+template <> //Work here
+template <>
+void particles_inner_computer_sedimenting_rods<double, long long>::compute_interaction_with_extra<6,6,9>(
+        const long long nb_particles,
+        const double pos_part[],
+        double rhs_part[],
+        const double rhs_part_extra[]) const{
+    //DEBUG_MSG("hello from compute_interaction_with_extra 9\n");
+    // call plain compute_interaction first
+    compute_interaction<6, 6>(nb_particles, pos_part, rhs_part);
+    const double ll2 = lambda*lambda;
+
+    // now add vorticity term
+    #pragma omp parallel for
+    for(long long idx_part = 0 ; idx_part < nb_particles ; ++idx_part){
+        const long long idx_part6 = idx_part*6 + 3;
+        const long long idx_part9 = idx_part*9;
+        rhs_part[idx_part6+IDXC_X] += (
+                pos_part[idx_part6+IDXC_Z]*(ll2*rhs_part_extra[idx_part9 + IDXC_DZ_X]-rhs_part_extra[idx_part9 + IDXC_DX_Z])
+              + pos_part[idx_part6+IDXC_Y]*(ll2*rhs_part_extra[idx_part9 + IDXC_DY_X]-rhs_part_extra[idx_part9 + IDXC_DX_Y])
+              + pos_part[idx_part6+IDXC_X]*(ll2-1)*rhs_part_extra[idx_part9 + IDXC_DX_X]) / (ll2+1);
+        rhs_part[idx_part6+IDXC_Y] += (
+                pos_part[idx_part6+IDXC_X]*(ll2*rhs_part_extra[idx_part9 + IDXC_DX_Y]-rhs_part_extra[idx_part9 + IDXC_DY_X])
+              + pos_part[idx_part6+IDXC_Z]*(ll2*rhs_part_extra[idx_part9 + IDXC_DZ_Y]-rhs_part_extra[idx_part9 + IDXC_DY_Z])
+              + pos_part[idx_part6+IDXC_Y]*(ll2-1)*rhs_part_extra[idx_part9 + IDXC_DY_Y]) / (ll2+1);
+        rhs_part[idx_part6+IDXC_Z] += (
+                pos_part[idx_part6+IDXC_Y]*(ll2*rhs_part_extra[idx_part9 + IDXC_DY_Z]-rhs_part_extra[idx_part9 + IDXC_DZ_Y])
+              + pos_part[idx_part6+IDXC_X]*(ll2*rhs_part_extra[idx_part9 + IDXC_DX_Z]-rhs_part_extra[idx_part9 + IDXC_DZ_X])
+              + pos_part[idx_part6+IDXC_Z]*(ll2-1)*rhs_part_extra[idx_part9 + IDXC_DZ_Z]) / (ll2+1);
+    }
+}
+
+
+// meant to be called AFTER executing the time-stepping operation.
+// once the particles have been moved, ensure that the orientation is a unit vector.
+template <>
+template <>
+void particles_inner_computer_sedimenting_rods<double, long long>::enforce_unit_orientation<6>(
+        const long long nb_particles,
+        double pos_part[]) const{
+
+    //DEBUG_MSG("hello from enforce_unit_orientation\n");
+    #pragma omp parallel for
+    for(long long idx_part = 0 ; idx_part < nb_particles ; ++idx_part){
+        const long long idx0 = idx_part*6 + 3;
+        // compute orientation size:
+        const double orientation_size = sqrt(
+                pos_part[idx0+IDXC_X]*pos_part[idx0+IDXC_X] +
+                pos_part[idx0+IDXC_Y]*pos_part[idx0+IDXC_Y] +
+                pos_part[idx0+IDXC_Z]*pos_part[idx0+IDXC_Z]);
+        // now renormalize
+        pos_part[idx0 + IDXC_X] /= orientation_size;
+        pos_part[idx0 + IDXC_Y] /= orientation_size;
+        pos_part[idx0 + IDXC_Z] /= orientation_size;
+    }
+}
+
+template
+void particles_inner_computer_sedimenting_rods<double, long long>::compute_interaction<6, 6>(
+        const long long nb_particles,
+        const double pos_part[],
+        double rhs_part[]) const;
\ No newline at end of file

projects/sedimenting_rods/particles_inner_computer_sedimenting_rods.hpp

+/******************************************************************************
+*                                                                             *
+*  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                                         *
+*                                                                             *
+******************************************************************************/
+
+
+
+#ifndef PARTICLES_INNER_COMPUTER_SEDIMENTING_RODS_HPP
+#define PARTICLES_INNER_COMPUTER_SEDIMENTING_RODS_HPP
+
+#include <cstring>
+#include <cassert>
+#include <iostream>
+
+template <class real_number, class partsize_t>
+class particles_inner_computer_sedimenting_rods{
+    bool isActive;
+    const real_number D_rho, u_sed;
+    const real_number beta0, beta1;
+    const real_number lambda;
+    const real_number lambda1;
+    const real_number lambda2;
+    const real_number lambda3;
+
+public:
+    explicit particles_inner_computer_sedimenting_rods(const real_number inLambda, const real_number inD_rho, const real_number inu_sed, const real_number in_beta0, const real_number in_beta1):
+        isActive(true),
+        D_rho(inD_rho),
+        u_sed(inu_sed),
+        beta0(in_beta0),
+        beta1(in_beta1),
+        lambda(inLambda),
+        lambda1(0),
+        lambda2(0),
+        lambda3(0)
+    {}
+
+    template <int size_particle_positions, int size_particle_rhs>
+    void compute_interaction(
+            const partsize_t nb_particles,
+            const real_number pos_part[],
+            real_number rhs_part[]) const;
+    // for given orientation and right-hand-side, recompute right-hand-side such
+    // that it is perpendicular to the current orientation.
+    // this is the job of the Lagrange multiplier terms, hence the
+    // "add_Lagrange_multipliers" name of the method.
+    template <int size_particle_positions, int size_particle_rhs>
+    void add_Lagrange_multipliers(
+            const partsize_t nb_particles,
+            const real_number pos_part[],
+            real_number rhs_part[]) const;
+    template <int size_particle_positions, int size_particle_rhs, int size_particle_rhs_extra>
+    void compute_interaction_with_extra(
+            const partsize_t nb_particles,
+            const real_number pos_part[],
+            real_number rhs_part[],
+            const real_number rhs_part_extra[]) const;
+    // meant to be called AFTER executing the time-stepping operation.
+    // once the particles have been moved, ensure that the orientation is a unit vector.
+    template <int size_particle_positions>
+    void enforce_unit_orientation(
+            const partsize_t nb_particles,
+            real_number pos_part[]) const;
+
+    bool isEnable() const {
+        return isActive;
+    }
+
+    void setEnable(const bool inIsActive) {
+        isActive = inIsActive;
+    }
+
+    double get_sedimenting_velocity () {
+        return u_sed;
+    }
+
+    void set_sedimenting_velocity (const double bla) {
+        this->u_sed = bla;
+    }
+
+    void set_D_rho (const double DRHO) {
+        this->D_rho = DRHO;
+    }
+
+    void set_beta0 (const double bla) {
+        this->beta0 = bla;
+    }
+
+    void set_beta1 (const double bla) {
+        this->beta1 = bla;
+    }
+
+    double get_D_rho () {
+        return D_rho;
+    }
+
+    double get_beta0 () {
+        return beta0;
+    }
+    double get_beta1 () {
+        return beta1;
+    }
+};
+
+#endif
\ No newline at end of file