diff --git a/projects/sedimenting_rods/GCDNS.py b/projects/sedimenting_rods/GCDNS.py
new file mode 100644
index 0000000000000000000000000000000000000000..670c910d270d3b5fb60aaf255b1b8f438cd0c77a
--- /dev/null
+++ b/projects/sedimenting_rods/GCDNS.py
@@ -0,0 +1,228 @@
+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()
diff --git a/projects/sedimenting_rods/README.rst b/projects/sedimenting_rods/README.rst
new file mode 100644
index 0000000000000000000000000000000000000000..b5dc14efcf3835a7cfe04e59c1a60adc05d04dbd
--- /dev/null
+++ b/projects/sedimenting_rods/README.rst
@@ -0,0 +1,9 @@
+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.
diff --git a/projects/sedimenting_rods/many_NSVEcylinder_collisions.cpp b/projects/sedimenting_rods/many_NSVEcylinder_collisions.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c231f2d3bf61b3ad8533f691025223551b65536d
--- /dev/null
+++ b/projects/sedimenting_rods/many_NSVEcylinder_collisions.cpp
@@ -0,0 +1,663 @@
+/**********************************************************************
+* *
+* 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 *
+* *
+**********************************************************************/
+
+#include <string>
+#include <cmath>
+#include "scope_timer.hpp"
+#include "particles/particles_sampling.hpp"
+#include "particles/p2p/p2p_computer.hpp"
+#include "particles/inner/particles_inner_computer.hpp"
+#include "particles/interpolation/field_tinterpolator.hpp"
+#include "particles/interpolation/particle_set.hpp"
+#include "particles/p2p/p2p_ghost_collisions.hpp"
+
+template <typename rnumber>
+int many_NSVEcylinder_collisions<rnumber>::initialize(void)
+{
+ TIMEZONE("many_NSVEcylinder_collisions::initialize");
+ this->NSVE<rnumber>::initialize();
+ this->many_ps.resize(this->many_integration_steps.size());
+ this->particles_output_writer_mpi.resize(this->many_integration_steps.size());
+ this->particles_sample_writer_mpi.resize(this->many_integration_steps.size());
+
+ // initialize particles
+ //force inner term to true.
+ //not sure why, since it's set to true in the constructor prototype
+ this->enable_inner = true;
+
+ for (unsigned int index = 0; index < this->many_ps.size(); index++)
+ {
+ DEBUG_MSG("in many_ps loop %d\n", index);
+ p2p_ghost_collisions<double, long long int> p2p_cc;
+ particles_inner_computer_sedimenting_rods<double, long long int> current_particles_inner_computer(this->tracers0_lambda[index], this->tracers0_D_rho[index], this->tracers0_u_sed[index], this->tracers0_beta0[index], this->tracers0_beta1[index]);
+ current_particles_inner_computer.setEnable(enable_inner);
+
+ DEBUG_MSG("Settling velocity is %g\n", current_particles_inner_computer.get_sedimenting_velocity() );
+ DEBUG_MSG("D_rho is %g\n", current_particles_inner_computer.get_D_rho() );
+ DEBUG_MSG("beta0 is %g\n", current_particles_inner_computer.get_beta0() );
+ DEBUG_MSG("beta1 is %g\n", current_particles_inner_computer.get_beta1() );
+ DEBUG_MSG("Is enable is %d\n", current_particles_inner_computer.isEnable() );
+
+ this->many_ps[index] = Template_double_for_if::evaluate<std::unique_ptr<abstract_particles_system_with_p2p<long long int, double, p2p_ghost_collisions<double, long long int>>>,
+ int, 1, 11, 1, // interpolation_size
+ int, 0, 3, 1, // spline_mode
+ particles_system_build_container<long long int, rnumber,FFTW,THREE,double,
+ p2p_ghost_collisions<double, long long int>,
+ particles_inner_computer_sedimenting_rods<double,long long int>,
+ 6,6>>(
+ this->tracers0_neighbours, // template iterator 1
+ this->tracers0_smoothness, // template iterator 2
+ this->fs->cvelocity,
+ this->fs->kk,
+ tracers0_integration_steps,
+ (long long int)nparticles[index],
+ this->get_current_fname(),
+ std::string("/tracers" + std::to_string(index) + "/state/") + std::to_string(this->iteration),
+ std::string("/tracers" + std::to_string(index) + "/rhs/") + std::to_string(this->iteration),
+ this->comm,
+ this->iteration+1,
+ std::move(p2p_cc),
+ std::move(current_particles_inner_computer),
+ this->tracers0_cutoff[index]);
+
+ //DEBUG_MSG("Settling velocity is %g\n", this->ps->getInnerComputer().get_sedimenting_velocity() );
+ //DEBUG_MSG("Is enable is %d\n", this->ps->getInnerComputer().isEnable() );
+
+ if (this->tracers0_lambda[index]>1.0){
+ this->many_ps[index]->getP2PComputer().setEnable(true);
+ this->many_ps[index]->getP2PComputer().set_cylinder();
+ this->many_ps[index]->getP2PComputer().set_cylinder_length(this->tracers0_cutoff[index]*(1.0-1.0/this->tracers0_lambda[index]));
+ this->many_ps[index]->getP2PComputer().set_cylinder_width(this->tracers0_cutoff[index] / this->tracers0_lambda[index]);
+ } else if (this->tracers0_lambda[index]==1.0){
+ this->many_ps[index]->getP2PComputer().set_sphere();
+ }
+
+ // initialize output objects
+ this->particles_output_writer_mpi[index] = new particles_output_hdf5<
+ long long int, double, 6>(
+ MPI_COMM_WORLD,
+ "tracers" + std::to_string(index),
+ nparticles[index],
+ this->many_integration_steps[index]);
+ this->particles_output_writer_mpi[index]->setParticleFileLayout(this->many_ps[index]->getParticleFileLayout());
+
+ }
+
+ /// allocate grad vel field
+ this->nabla_u = new field<rnumber, FFTW, THREExTHREE>(
+ this->nx, this->ny, this->nz,
+ this->comm,
+ this->fs->cvorticity->fftw_plan_rigor);
+ /// open collisions file
+ 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 cols_file cache I got %d\n", cache_err);
+ this->cols_file = H5Fopen(
+ (this->simname + "_collisions.h5").c_str(),
+ H5F_ACC_RDWR,
+ fapl);
+ this->cols_file_PP = H5Fopen(
+ (this->simname + "_collisions_PP.h5").c_str(),
+ H5F_ACC_RDWR,
+ fapl);
+ this->collision_pairs_old.resize(this->many_ps.size());
+ this->collision_pairs_new.resize(this->many_ps.size());
+ this->collision_pairs_same.resize(this->many_ps.size());
+ this->collision_pairs_time_0.resize(this->many_ps.size());
+ this->collision_pairs_time_1.resize(this->many_ps.size());
+ for (unsigned int index = 0; index < this->many_ps.size(); index++)
+ {
+ DEBUG_MSG("index %d\n", index);
+ reset_collision_lists(index);
+ }
+ }
+
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int many_NSVEcylinder_collisions<rnumber>::step(void)
+{
+ TIMEZONE("many_NSVEcylinder_collisions::step");
+ this->fs->compute_velocity(this->fs->cvorticity);
+
+ compute_gradient(
+ this->fs->kk,
+ this->fs->cvelocity,
+ this->nabla_u);
+ this->nabla_u->ift();
+ this->fs->cvelocity->ift(); // needed before completeloop
+ //std::unique_ptr<double[]> sampled_vorticity(new double[9*this->ps->getLocalNbParticles()]);
+ //std::fill_n(sampled_vorticity.get(), 9*this->ps->getLocalNbParticles(), 0);
+ //this->ps->sample_compute_field(*this->nabla_u, sampled_vorticity.get());
+ //*this->tmp_vec_field = this->fs->cvorticity->get_cdata();
+ //this->tmp_vec_field->ift();
+
+ for (unsigned int index = 0; index < this->many_ps.size(); index++)
+ {
+ DEBUG_MSG(("step for tracers"+std::to_string(index)+"\n").c_str());
+ this->many_ps[index]->getP2PComputer().reset_collision_pairs();
+ DEBUG_MSG(("step for tracers"+std::to_string(index)+", after reset_collision_pairs\n").c_str());
+ this->many_ps[index]->completeLoopWithExtraField(this->dt, *this->nabla_u);
+ DEBUG_MSG(("step for tracers"+std::to_string(index)+", after completeLoopWithExtraField\n").c_str());
+ }
+
+ this->NSVE<rnumber>::step();
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int many_NSVEcylinder_collisions<rnumber>::write_checkpoint(void)
+{
+ TIMEZONE("many_NSVEcylinder_collisions::write_checkpoint");
+ this->NSVE<rnumber>::write_checkpoint();
+ for (unsigned int index = 0; index < this->many_ps.size(); index++)
+ {
+ this->particles_output_writer_mpi[index]->open_file(this->fs->get_current_fname());
+ // TODO P2P write particle data too
+ this->particles_output_writer_mpi[index]->template save<6>(
+ this->many_ps[index]->getParticlesState(),
+ this->many_ps[index]->getParticlesRhs(),
+ this->many_ps[index]->getParticlesIndexes(),
+ this->many_ps[index]->getLocalNbParticles(),
+ this->fs->iteration);
+ this->particles_output_writer_mpi[index]->close_file();
+ }
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int many_NSVEcylinder_collisions<rnumber>::finalize(void)
+{
+ TIMEZONE("many_NSVEcylinder_collisions::finalize");
+ if (this->myrank == 0){
+ H5Fclose(this->cols_file);
+ H5Fclose(this->cols_file_PP);
+ }
+ delete this->nabla_u;
+ for (unsigned int index = 0; index < this->many_ps.size(); index++)
+ {
+ delete this->particles_output_writer_mpi[index];
+ }
+ this->NSVE<rnumber>::finalize();
+
+ return EXIT_SUCCESS;
+}
+
+/** \brief Compute fluid stats and sample particle data.
+ */
+
+template <typename rnumber>
+int many_NSVEcylinder_collisions<rnumber>::do_stats()
+{
+ bool save_pair;
+ long long int indx_0_0, indx_0_1, indx_1_0, indx_1_1;
+ TIMEZONE("many_NSVEcylinder_collisions::do_stats");
+ /// perform fluid stats
+ this->NSVE<rnumber>::do_stats();
+
+ // compute velocity, copy to tmp field
+ //this->fs->compute_velocity(this->fs->cvorticity);
+ //*(this->tmp_vec_field) = *this->fs->cvelocity;
+
+ /// sample velocity gradient
+ /// fs->cvelocity should contain the velocity in Fourier space
+ this->fs->compute_velocity(this->fs->cvorticity);
+ compute_gradient(
+ this->fs->kk,
+ this->fs->cvelocity,
+ this->nabla_u);
+ this->nabla_u->ift();
+ // Dont need acceleration in the moment //
+ ///// compute acceleration and sample it
+ //this->fs->compute_Lagrangian_acceleration(this->tmp_vec_field);
+ // move velocity to real space
+ this->fs->cvelocity->ift();
+
+ // Carry out collision measurement
+ if (this->enable_p2p) {
+ for (unsigned int index = 0; index < this->many_ps.size(); index++)
+ {
+ /* DEBUG_MSG("Saving collisions at index %d\n",this->many_ps[index]->get_step_idx()-1); */
+ /* Here the collisions are saved both at the end and at the beginning of a subsequent simulation. This produces an error message. */
+ DEBUG_MSG(("merge for tracers"+std::to_string(index)+"\n").c_str());
+ this->many_ps[index]->getP2PComputer().MPI_merge(
+ this->comm,
+ this->myrank,
+ this->nprocs);
+
+ if (this->myrank == 0){
+ /* Reset the collision lists old, same, and new */
+ this->reset_collision_lists(index);
+ /* Copy the collision list from the current time step onto collision_pairs_time_1 */
+ this->collision_pairs_time_1[index] = this->many_ps[index]->getP2PComputer().get_collision_pairs(
+ this->comm,
+ this->myrank,
+ this->nprocs);
+ /* Get the collisions that dissapeared in this time step and those that remain the same */
+ for (unsigned int pair_index_0 = 0; pair_index_0 < this->collision_pairs_time_0[index].size()/2; pair_index_0++)
+ {
+ indx_0_0 = this->collision_pairs_time_0[index][pair_index_0*2];
+ indx_0_1 = this->collision_pairs_time_0[index][pair_index_0*2+1];
+ save_pair = true;
+ for (unsigned int pair_index_1 = 0; pair_index_1 < this->collision_pairs_time_1[index].size()/2; pair_index_1++)
+ {
+ indx_1_0 = this->collision_pairs_time_1[index][pair_index_1*2];
+ indx_1_1 = this->collision_pairs_time_1[index][pair_index_1*2+1];
+ if (indx_0_0==indx_1_0 and indx_0_1==indx_1_1)
+ {
+ save_pair = false;
+ this->collision_pairs_same[index].push_back(indx_0_0);
+ this->collision_pairs_same[index].push_back(indx_0_1);
+ break;
+ }
+ }
+ if (save_pair){
+ this->collision_pairs_old[index].push_back(indx_0_0);
+ this->collision_pairs_old[index].push_back(indx_0_1);
+ }
+ }
+ /* Get the new collisions from this time step */
+ for (unsigned int pair_index_0 = 0; pair_index_0 < this->collision_pairs_time_1[index].size()/2; pair_index_0++)
+ {
+ indx_0_0 = this->collision_pairs_time_1[index][pair_index_0*2];
+ indx_0_1 = this->collision_pairs_time_1[index][pair_index_0*2+1];
+ save_pair = true;
+ for (unsigned int pair_index_1 = 0; pair_index_1 < this->collision_pairs_time_0[index].size()/2; pair_index_1++)
+ {
+ indx_1_0 = this->collision_pairs_time_0[index][pair_index_1*2];
+ indx_1_1 = this->collision_pairs_time_0[index][pair_index_1*2+1];
+ if (indx_0_0==indx_1_0 and indx_0_1==indx_1_1)
+ {
+ save_pair = false;
+ break;
+ }
+ }
+ if (save_pair){
+ this->collision_pairs_new[index].push_back(indx_0_0);
+ this->collision_pairs_new[index].push_back(indx_0_1);
+ }
+ }
+ DEBUG_MSG(("save cols for tracers"+std::to_string(index)+"\n").c_str());
+ // this->cols_file is only defined for rank 0,
+ // so only rank 0 can write the pairs.
+ // this is also a reason for calling MPI_merge beforehand.
+ std::string dset_name = (
+ "/collisions"
+ + std::to_string(index)
+ + "/pair_list_"
+ + std::to_string(this->many_ps[index]->get_step_idx()-1));
+ // if the dataset already exists, it's because we're doing the statistics for
+ // first iteration in current run. this is fine for spectra etc (what happens for example
+ // in NSVE::do_stats), where the same element of some big array gets rewritten.
+ // However, it's not fine for "write_particle_ID_pairs", because this method creates
+ // a new dataset. Which raises an error when the dataset already exists.
+ // Hence I chose to write the current data with "extra", for redundancy checks.
+ // Although I assume the "extra" dataset will be empty since timestep has not run
+ // and no collisions have been counted.
+ if (H5Lexists(this->cols_file, dset_name.c_str(), H5P_DEFAULT))
+ dset_name += "_extra";
+ hdf5_tools::write_particle_ID_pairs_with_single_rank<long long>(
+ this->collision_pairs_time_1[index],
+ this->cols_file,
+ dset_name.c_str());
+ /* Save the new, old, and same collisions */
+ DEBUG_MSG(("save new cols for tracers"+std::to_string(index)+"\n").c_str());
+ dset_name = (
+ "/new_collisions"
+ + std::to_string(index)
+ + "/pair_list_"
+ + std::to_string(this->many_ps[index]->get_step_idx()-1));
+ if (H5Lexists(this->cols_file, dset_name.c_str(), H5P_DEFAULT))
+ dset_name += "_extra";
+ hdf5_tools::write_particle_ID_pairs_with_single_rank<long long>(
+ this->collision_pairs_new[index],
+ this->cols_file,
+ dset_name.c_str());
+ hdf5_tools::write_particle_ID_pairs_with_single_rank<long long>(
+ this->collision_pairs_new[index],
+ this->cols_file_PP,
+ dset_name.c_str());
+ DEBUG_MSG(("save old cols for tracers"+std::to_string(index)+"\n").c_str());
+ dset_name = (
+ "/old_collisions"
+ + std::to_string(index)
+ + "/pair_list_"
+ + std::to_string(this->many_ps[index]->get_step_idx()-1));
+ if (H5Lexists(this->cols_file, dset_name.c_str(), H5P_DEFAULT))
+ dset_name += "_extra";
+ hdf5_tools::write_particle_ID_pairs_with_single_rank<long long>(
+ this->collision_pairs_old[index],
+ this->cols_file,
+ dset_name.c_str());
+ DEBUG_MSG(("save same cols for tracers"+std::to_string(index)+"\n").c_str());
+ dset_name = (
+ "/same_collisions"
+ + std::to_string(index)
+ + "/pair_list_"
+ + std::to_string(this->many_ps[index]->get_step_idx()-1));
+ if (H5Lexists(this->cols_file, dset_name.c_str(), H5P_DEFAULT))
+ dset_name += "_extra";
+ hdf5_tools::write_particle_ID_pairs_with_single_rank<long long>(
+ this->collision_pairs_same[index],
+ this->cols_file,
+ dset_name.c_str());
+ DEBUG_MSG(("save stats for tracers"+std::to_string(index)+"\n").c_str());
+ /* Copy the collision list from the current time step onto collision_pairs_time_0 */
+ this->collision_pairs_time_0[index]=this->collision_pairs_time_1[index];
+ DEBUG_MSG(("done saving stats for tracers"+std::to_string(index)+"\n").c_str());
+ }
+
+ {
+ // this block for subsampling
+ // declare temporary particle set for subset sampling
+ DEBUG_MSG_WAIT(this->comm, "00 hello from subsample block\n");
+ particle_set<6, 3, 2> psubset(
+ this->fs->cvelocity->rlayout,
+ this->dkx,
+ this->dky,
+ this->dkz,
+ this->tracers0_cutoff[index]);
+ DEBUG_MSG_WAIT(this->comm, "01 hello from subsample block\n");
+
+ // first extract unique IDs of colliding particles
+ int nparticles_to_sample;
+ std::vector<long long int> indices_to_sample;
+ if (this->myrank == 0)
+ {
+ indices_to_sample.resize(this->collision_pairs_new[index].size());
+ // direct copy
+ for (unsigned ii = 0; ii < indices_to_sample.size(); ii++)
+ indices_to_sample[ii] = this->collision_pairs_new[index][ii];
+ // sort in place
+ std::sort(indices_to_sample.begin(), indices_to_sample.end());
+ // remove adjacent duplicates
+ auto last = std::unique(indices_to_sample.begin(), indices_to_sample.end());
+ // clean up undefined values
+ indices_to_sample.erase(last, indices_to_sample.end());
+ // now indices_to_sample contains unique particle labels.
+ nparticles_to_sample = indices_to_sample.size();
+ }
+ DEBUG_MSG_WAIT(this->comm, "02 hello from subsample block\n");
+
+ // tell all processes how many particles there are
+ AssertMpi(MPI_Bcast(
+ &nparticles_to_sample,
+ 1,
+ MPI_INT,
+ 0,
+ this->comm));
+ DEBUG_MSG_WAIT(this->comm, "03 hello from subsample block\n");
+ if (nparticles_to_sample > 0)
+ {
+ DEBUG_MSG_WAIT(this->comm, "hello, there are %d particles to sample\n", nparticles_to_sample);
+ // tell all processes which particles are in the subset
+ if (this->myrank != 0)
+ indices_to_sample.resize(nparticles_to_sample);
+ AssertMpi(MPI_Bcast(
+ &indices_to_sample.front(),
+ nparticles_to_sample,
+ MPI_LONG,
+ 0,
+ this->comm));
+ // now initialize particle set
+ psubset.init_as_subset_of(
+ *this->many_ps[index],
+ indices_to_sample);
+ DEBUG_MSG_WAIT(this->comm, "hello, psubset has %ld particles out of %ld\n",
+ psubset.getLocalNumberOfParticles(),
+ psubset.getTotalNumberOfParticles());
+ // write labels of colliding particles
+ // open subset file only with rank 0
+ // START WARNING: subset file opens "colliding_particle_subset"
+ std::string file_name = this->simname + "_colliding_particle_subset.h5";
+ if (this->myrank == 0)
+ this->subset_file = H5Fopen(
+ file_name.c_str(),
+ H5F_ACC_RDWR,
+ H5P_DEFAULT);
+ else
+ this->subset_file = -1;
+ DEBUG_MSG_WAIT(this->comm, "psubset check 00\n");
+ hdf5_tools::gather_and_write_with_single_rank(
+ this->myrank,
+ 0,
+ this->comm,
+ psubset.getParticleIndices(),
+ psubset.getLocalNumberOfParticles(),
+ this->subset_file,
+ "subset_tracers"+ std::to_string(index) + "/labels/" + std::to_string(this->iteration));
+ // write positions of colliding particles
+ DEBUG_MSG_WAIT(this->comm, "psubset check 01\n");
+ std::unique_ptr<double[]> position = psubset.extractFromParticleState(0, 3);
+ hdf5_tools::gather_and_write_with_single_rank(
+ this->myrank,
+ 0,
+ this->comm,
+ position.get(),
+ psubset.getLocalNumberOfParticles()*3,
+ this->subset_file,
+ "subset_tracers"+ std::to_string(index) + "/position/" + std::to_string(this->iteration));
+ delete[] position.release();
+ // write orientations of colliding particles
+ DEBUG_MSG_WAIT(this->comm, "psubset check 02\n");
+ std::unique_ptr<double[]> orientation = psubset.extractFromParticleState(3, 6);
+ hdf5_tools::gather_and_write_with_single_rank(
+ this->myrank,
+ 0,
+ this->comm,
+ orientation.get(),
+ psubset.getLocalNumberOfParticles()*3,
+ this->subset_file,
+ "subset_tracers"+ std::to_string(index) + "/orientation/" + std::to_string(this->iteration));
+ if (this->myrank == 0)
+ H5Fclose(this->subset_file);
+ DEBUG_MSG_WAIT(this->comm, "psubset check 03\n");
+ MPI_Barrier(this->comm);
+ // END WARNING: subset file "colliding_particle_subset"
+ delete[] orientation.release();
+ // START WARNING these calls use "colliding_particle_subset"
+ // now write velocity sample
+ psubset.writeSample(
+ this->fs->cvelocity,
+ this->simname + "_colliding_particle_subset.h5",
+ "subset_tracers"+ std::to_string(index),
+ "velocity",
+ this->iteration);
+ // now write gradient sample
+ psubset.writeSample(
+ this->nabla_u,
+ this->simname + "_colliding_particle_subset.h5",
+ "subset_tracers"+ std::to_string(index),
+ "velocity_gradient",
+ this->iteration);
+ // END WARNING
+ }
+ DEBUG_MSG_WAIT(this->comm, "04 hello from subsample block\n");
+ }
+ if (this->myrank == 0)
+ {
+ /* Reset the collision lists old, same, and new */
+ this->reset_collision_lists(index);
+ }
+ }
+ }
+
+ /// check if particle stats should be performed now;
+ /// if not, exit method.
+ if (!(this->iteration % this->niter_part == 0))
+ {
+ return EXIT_SUCCESS;
+ }
+
+
+ DEBUG_MSG_WAIT(this->comm, "before sample write\n");
+ for (unsigned int index = 0; index < this->many_ps.size(); index++)
+ {
+ this->particles_sample_writer_mpi[index] = new particles_output_sampling_hdf5<
+ long long int, double, double, 3>(
+ MPI_COMM_WORLD,
+ this->many_ps[index]->getGlobalNbParticles(),
+ (this->simname + "_particles.h5"),
+ "tracers" + std::to_string(index),
+ "position/0");
+
+ this->particles_sample_writer_mpi[index]->setParticleFileLayout(this->many_ps[index]->getParticleFileLayout());
+ DEBUG_MSG(("sampling tracers"+std::to_string(index)+"\n").c_str());
+ /// allocate temporary data array
+ /// initialize pdata0 with the positions, and pdata1 with the orientations
+ std::unique_ptr<double[]> pdata0 = this->many_ps[index]->extractParticlesState(0, 3);
+ std::unique_ptr<double[]> pdata1 = this->many_ps[index]->extractParticlesState(3, 6);
+ std::unique_ptr<double[]> pdata2(new double[9*this->many_ps[index]->getLocalNbParticles()]);
+ DEBUG_MSG(("sampling position for tracers"+std::to_string(index)+"\n").c_str());
+ /// sample position
+ this->particles_sample_writer_mpi[index]->template save_dataset<3>(
+ "tracers" + std::to_string(index),
+ "position",
+ pdata0.get(), // we need to use pdata0.get() here, because it's 3D, whereas getParticlesState may give something else
+ &pdata0,
+ this->many_ps[index]->getParticlesIndexes(),
+ this->many_ps[index]->getLocalNbParticles(),
+ this->many_ps[index]->get_step_idx()-1);
+ DEBUG_MSG(("sampling orientation for tracers"+std::to_string(index)+"\n").c_str());
+ /// sample orientation
+ this->particles_sample_writer_mpi[index]->template save_dataset<3>(
+ "tracers" + std::to_string(index),
+ "orientation",
+ pdata0.get(),
+ &pdata1,
+ this->many_ps[index]->getParticlesIndexes(),
+ this->many_ps[index]->getLocalNbParticles(),
+ this->many_ps[index]->get_step_idx()-1);
+ DEBUG_MSG(("sampling velocity for tracers"+std::to_string(index)+"\n").c_str());
+ /// sample velocity
+ /// from now on, we need to clean up data arrays before interpolation
+ std::fill_n(pdata1.get(), 3*this->many_ps[index]->getLocalNbParticles(), 0);
+ DEBUG_MSG(("filled pdata1 for tracers"+std::to_string(index)+"\n").c_str());
+ this->many_ps[index]->sample_compute_field(*this->fs->cvelocity, pdata1.get());
+ DEBUG_MSG(("sampled for tracers"+std::to_string(index)+"\n").c_str());
+ this->particles_sample_writer_mpi[index]->template save_dataset<3>(
+ "tracers" + std::to_string(index),
+ "velocity",
+ pdata0.get(),
+ &pdata1,
+ this->many_ps[index]->getParticlesIndexes(),
+ this->many_ps[index]->getLocalNbParticles(),
+ this->many_ps[index]->get_step_idx()-1);
+ DEBUG_MSG(("sampling velocity gradient for tracers"+std::to_string(index)+"\n").c_str());
+ std::fill_n(pdata2.get(), 9*this->many_ps[index]->getLocalNbParticles(), 0);
+ this->many_ps[index]->sample_compute_field(*this->nabla_u, pdata2.get());
+ this->particles_sample_writer_mpi[index]->template save_dataset<9>(
+ "tracers" + std::to_string(index),
+ "velocity_gradient",
+ pdata0.get(),
+ &pdata2,
+ this->many_ps[index]->getParticlesIndexes(),
+ this->many_ps[index]->getLocalNbParticles(),
+ this->many_ps[index]->get_step_idx()-1);
+
+ //Dont need accelaration at the moment.
+// std::fill_n(pdata1.get(), 3*this->many_ps[index]->getLocalNbParticles(), 0);
+// this->many_ps[index]->sample_compute_field(*this->tmp_vec_field, pdata1.get());
+// this->particles_sample_writer_mpi[index]->template save_dataset<3>(
+// "tracers" + std::to_string(index),
+// "acceleration",
+// pdata0.get(),
+// &pdata1,
+// this->many_ps[index]->getParticlesIndexes(),
+// this->many_ps[index]->getLocalNbParticles(),
+// this->many_ps[index]->get_step_idx()-1);
+ delete this->particles_sample_writer_mpi[index];
+ // deallocate temporary data array
+ delete[] pdata0.release();
+ delete[] pdata1.release();
+ delete[] pdata2.release();
+ }
+ DEBUG_MSG_WAIT(this->comm, "after sample write\n");
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int many_NSVEcylinder_collisions<rnumber>::read_parameters(void)
+{
+ TIMEZONE("many_NSVEcylinder_collisions::read_parameters");
+ this->NSVE<rnumber>::read_parameters();
+ hid_t parameter_file = H5Fopen((this->simname + ".h5").c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+ this->niter_part = hdf5_tools::read_value<int>(parameter_file, "parameters/niter_part");
+ DEBUG_MSG("Loading nparticles");
+ this->nparticles = hdf5_tools::read_vector<int>(
+ parameter_file,
+ "parameters/nparticles");
+ //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");
+ this->enable_p2p = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_enable_p2p");
+ this->enable_inner = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_enable_inner");
+ int tval = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_enable_vorticity_omega");
+ this->enable_vorticity_omega = tval;
+ DEBUG_MSG("tracers0_enable_vorticity_omega = %d, this->enable_vorticity_omega = %d\n",
+ tval, this->enable_vorticity_omega);
+ DEBUG_MSG("Loading tracers0_cutoff \n");
+ this->tracers0_cutoff = hdf5_tools::read_vector<double>(
+ parameter_file,
+ "parameters/tracers0_cutoff");
+ DEBUG_MSG("Loading tracers0_lambda \n");
+ this->tracers0_lambda = hdf5_tools::read_vector<double>(
+ parameter_file,
+ "parameters/tracers0_lambda");
+ DEBUG_MSG("Loading tracers0_D_rho \n");
+ this->tracers0_D_rho = hdf5_tools::read_vector<double>(
+ parameter_file,
+ "parameters/tracers0_D_rho");
+ DEBUG_MSG("Loading many_integration_steps \n");
+ this->many_integration_steps = hdf5_tools::read_vector<int>(
+ parameter_file,
+ "parameters/many_integration_steps");
+ DEBUG_MSG("Loading tracers0_u_sed \n");
+ this->tracers0_u_sed = hdf5_tools::read_vector<double>(
+ parameter_file,
+ "parameters/tracers0_u_sed");
+ DEBUG_MSG("Loading tracers0_beta0 \n");
+ this->tracers0_beta0 = hdf5_tools::read_vector<double>(
+ parameter_file,
+ "parameters/tracers0_beta0");
+ DEBUG_MSG("Loading tracers0_beta1 \n");
+ this->tracers0_beta1 = hdf5_tools::read_vector<double>(
+ parameter_file,
+ "parameters/tracers0_beta1");
+ 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");
+ H5Fclose(parameter_file);
+ MPI_Barrier(this->comm);
+ assert(this->many_integration_steps.size() == this->tracers0_D_rho.size());
+ return EXIT_SUCCESS;
+}
+
+
+template class many_NSVEcylinder_collisions<float>;
+template class many_NSVEcylinder_collisions<double>;
diff --git a/projects/sedimenting_rods/many_NSVEcylinder_collisions.hpp b/projects/sedimenting_rods/many_NSVEcylinder_collisions.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d89526c969c931a61a204d6365e406bee93fc83a
--- /dev/null
+++ b/projects/sedimenting_rods/many_NSVEcylinder_collisions.hpp
@@ -0,0 +1,122 @@
+/**********************************************************************
+* *
+* 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
diff --git a/projects/sedimenting_rods/particles_inner_computer_sedimenting_rods.cpp b/projects/sedimenting_rods/particles_inner_computer_sedimenting_rods.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..101eaf3631550898d4a97bfa1380205e25a76e08
--- /dev/null
+++ b/projects/sedimenting_rods/particles_inner_computer_sedimenting_rods.cpp
@@ -0,0 +1,239 @@
+/******************************************************************************
+* *
+* 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
diff --git a/projects/sedimenting_rods/particles_inner_computer_sedimenting_rods.hpp b/projects/sedimenting_rods/particles_inner_computer_sedimenting_rods.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..6fa808ff7cc0c663088a66cc2cd8f8cff0ac544d
--- /dev/null
+++ b/projects/sedimenting_rods/particles_inner_computer_sedimenting_rods.hpp
@@ -0,0 +1,123 @@
+/******************************************************************************
+* *
+* 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