diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000..5172429 Binary files /dev/null and b/.DS_Store differ diff --git a/Makefile b/Makefile index a67c606..c9f4488 100644 --- a/Makefile +++ b/Makefile @@ -2,9 +2,9 @@ # Path to the GSL, Boost C++ and CGAL libraries - must be set by user (only if they aren't installed in the default system path) -- (NOTE: You must add only the directory where the libraries are installed, the program will add the '/lib' and '/include' parts automatically); C++ compiler - preferably a version that supports OpenMP -GSL_PATH = /net/plato/data/users/cautun/Programs/stow/ -BOOST_PATH = -CGAL_PATH = +GSL_PATH = /cosma/local/gsl/2.4 +BOOST_PATH = /cosma/local/boost/gnu_7.3.0/1_67_0 +CGAL_PATH = /cosma/home/dphlss/cautun/Programs/stow CC = g++ # set the following if you have installed the HDF5 library and would like to read in HDF5 gadget files HDF5_PATH = @@ -105,9 +105,7 @@ endif COMPILE_FLAGS = -frounding-math -O3 -fopenmp -DNDEBUG $(OPTIONS) DTFE_INC = $(INCLUDES) # the following libraries should work in most cases -DTFE_LIB = -rdynamic $(LIBRARIES) -lCGAL -lboost_thread -lboost_filesystem -lboost_program_options -lgsl -lgslcblas -lm -lboost_system -# the following libraries work on "Fedora 15" -# DTFE_LIB = -rdynamic $(LIBRARIES) -lCGAL -lboost_thread-mt -lboost_filesystem -lboost_program_options -lgsl -lgslcblas -lgmp -lboost_system +DTFE_LIB = -rdynamic $(LIBRARIES) -lCGAL -lboost_thread -lboost_filesystem -lboost_program_options -lgmp -lgsl -lgslcblas -lm -lboost_system diff --git a/src/CIC_interpolation.cc b/src/CIC_interpolation.cc index bd529ad..a70ccde 100644 --- a/src/CIC_interpolation.cc +++ b/src/CIC_interpolation.cc @@ -44,7 +44,7 @@ void CIC_interpolation_regular_grid(vector &particles, -/* This function interpolates the density and velocity to grid using the CIC (cloud in cell) method. */ +/* This function interpolates the density and velocity to grid using the CIC (triangular shape cloud) method. */ void CIC_interpolation(vector *particles, vector &samples, User_options &userOptions, @@ -79,47 +79,184 @@ void CIC_interpolation_regular_grid(vector &particles, // allocate memory for the results size_t const reserveSize = (NO_DIM==2) ? nGrid[0]*nGrid[1] : nGrid[0]*nGrid[1]*nGrid[2]; q->density.assign( reserveSize, Real(0.) ); - if ( userOptions.aField.velocity ) - q->velocity.assign( reserveSize, Pvector::zero() ); + q->velocity.assign( reserveSize, Pvector::zero() ); - // get the grid spacing - Box box = userOptions.region; - Real dx[NO_DIM]; + // get the grid spacing and find the extended box of particles that give contribution in the region of interest + Box box = userOptions.region, + outerBox = userOptions.region, + innerBox = userOptions.region; + Real dx[NO_DIM], outerPadding[2*NO_DIM], innerPadding[2*NO_DIM]; for (int i=0; i innerParticles, outerParticles; for (vectorIterator it=particles.begin(); it!=particles.end(); ++it) - { - int cell[NO_DIM]; - bool validCell = true; - for (int j=0; jposition(j)-box[2*j])/dx[j] ) ); - if ( cell[j]<0 or cell[j]>=nGrid[j] ) - validCell = false; - } - if ( not validCell ) continue; + if ( innerBox.isParticleInBox(*it) ) + innerParticles.push_back( it ); // keep track of the particles that give contributions only inside the region of interest + else if ( outerBox.isParticleInBox(*it) ) + outerParticles.push_back( it ); // these particles give contributions also outside the region of interest + + + + // loop over all the particles in the inner box + for (list< vectorIterator >::iterator it=innerParticles.begin(); it!=innerParticles.end(); ++it) + { + int cell[NO_DIM][3]; // the grid coordinates of the cell that get a contribution from the particle + Real weight[NO_DIM][3]; // the weight associated to the particle for each cell that it contributes to + for (int j=0; jposition(j) - box[2*j] ) / dx[j]; + cell[j][0] = int(floor( temp )) - 1; + cell[j][1] = cell[j][0] + 1; + cell[j][2] = cell[j][1] + 1; + + temp -= cell[j][1] + 0.5; //this gives the distance of the particle with respect to the center of the density grid in which the particle is located + + // compute the weight associated to this direction + if ( temp < 0.) + { + weight[j][0] = -temp; + weight[j][1] = 1 + temp; + weight[j][2] = 0.; + } + else + { + weight[j][0] = 0.; + weight[j][1] = 1 - temp; + weight[j][2] = temp; + } + + } + + // get the density contribution of the particle to the neighboring cells +#if NO_DIM==2 + for (int i1=0; i1<3; ++i1) + for (int i2=0; i2<3; ++i2) + { + int index = cell[0][i1] * nGrid[1] + cell[1][i2]; + Real result = (*it)->weight() * weight[0][i1] * weight[1][i2]; + q->density[index] += result; + q->velocity[index] += (*it)->velocity() * result; + } +#elif NO_DIM==3 + for (int i1=0; i1<3; ++i1) + for (int i2=0; i2<3; ++i2) + for (int i3=0; i3<3; ++i3) + { + int index = cell[0][i1] * nGrid[1]*nGrid[2] + cell[1][i2] * nGrid[2] + cell[2][i3]; + Real result = (*it)->weight() * weight[0][i1] * weight[1][i2] * weight[2][i3]; + q->density[index] += result; + q->velocity[index] += ( (*it)->velocity() * result ); + } +#endif + } + + + // now loop over the particles on the boundary - must check that all neighbors are valid cells + for (list< vectorIterator >::iterator it=outerParticles.begin(); it!=outerParticles.end(); ++it) + { + int cell[NO_DIM][3]; // the grid coordinates of the cell that get a contribution from the particle + int cellCount[NO_DIM]; // counts how many valid cells are along each dimension that get contribution from the particle + for (int i=0; iposition(j) - box[2*j] ) / dx[j]; + int tempInt = int(floor( temp )); + temp -= tempInt + 0.5; //this gives the distance of the particle with respect to the center of the density grid in which the particle is located + if ( temp<-0.5 or temp>0.5 ) cout << temp << "\t" << j << "\t" << tempInt << "\n"; + + // check that all neighbors must be valid cells + if ( tempInt==-1 ) + { + cell[j][0] = 0; cellCount[j] = 1; + weight[j][0] = temp; + } + else if ( tempInt==0 ) + { + cell[j][0] = 0; cell[j][1] = 1; cellCount[j] = 2; + weight[j][0] = 1 - temp; + weight[j][1] = temp; + } + else if ( tempInt==int(nGrid[j])-1 ) + { + cell[j][0] = nGrid[j]-2; cell[j][1] = nGrid[j]-1; cellCount[j] = 2; + weight[j][0] = - temp; + weight[j][1] = 1 + temp; + } + else if ( tempInt==int(nGrid[j]) ) + { + cell[j][0] = nGrid[j]-1; cellCount[j] = 1; + weight[j][0] = - temp; + } + else if ( not( tempInt<-1 or tempInt>int(nGrid[j]) ) ) + { + cell[j][0] = tempInt-1; + cell[j][1] = tempInt; + cell[j][2] = tempInt+1; + cellCount[j] = 3; + + if ( temp < 0.) + { + weight[j][0] = -temp; + weight[j][1] = 1 + temp; + weight[j][2] = 0.; + } + else + { + weight[j][0] = 0.; + weight[j][1] = 1 - temp; + weight[j][2] = temp; + } + } + } + + // get the density contribution of the particle to the neighboring cells #if NO_DIM==2 - int index = cell[0] * nGrid[1] + cell[1]; + for (int i1=0; i1weight() * weight[0][i1] * weight[1][i2]; + q->density[index] += result; + q->velocity[index] += (*it)->velocity() * result; + } #elif NO_DIM==3 - int index = cell[0] * nGrid[1]*nGrid[2] + cell[1] * nGrid[2] + cell[2]; + for (int i1=0; i1weight() * weight[0][i1] * weight[1][i2] * weight[2][i3]; + q->density[index] += result; + q->velocity[index] += (*it)->velocity() * result; + } #endif - q->density[index] += it->weight(); - if ( userOptions.aField.velocity ) - q->velocity[index] += it->velocity() * it->weight(); - } + } // divide the momentum by the mass in the cells if ( userOptions.aField.velocity ) + { for (size_t i=0; idensity[i]!=Real(0.) ) q->velocity[i] /= q->density[i]; else q->velocity[i] = Pvector::zero(); + } + else + q->velocity.clear(); // normalize the density to average background density if ( userOptions.aField.density ) @@ -136,42 +273,5 @@ void CIC_interpolation_regular_grid(vector &particles, } -/* This function counts how many particles are in each cell of a grid using the CIC method. */ -void CIC_particle_count(vector &particles, - size_t const *nGrid, - Box box, - vector *counts) -{ - // allocate memory for the results - size_t const reserveSize = (NO_DIM==2) ? nGrid[0]*nGrid[1] : nGrid[0]*nGrid[1]*nGrid[2]; - counts->assign( reserveSize, int(0) ); - - // get the grid spacing - Real dx[NO_DIM]; - for (int i=0; iposition(j)-box[2*j])/dx[j] ) ); - if ( cell[j]<0 or cell[j]>=nGrid[j] ) - validCell = false; - } - if ( not validCell ) continue; -#if NO_DIM==2 - int index = cell[0] * nGrid[1] + cell[1]; -#elif NO_DIM==3 - int index = cell[0] * nGrid[1]*nGrid[2] + cell[1] * nGrid[2] + cell[2]; -#endif - (*counts)[index] += 1; - } -} - diff --git a/src/DTFE.cpp b/src/DTFE.cpp index 5257777..c37ca8f 100644 --- a/src/DTFE.cpp +++ b/src/DTFE.cpp @@ -41,6 +41,7 @@ using namespace std; #include "user_options.cc" #include "quantities.cc" +#include "NGP_interpolation.cc" #include "CIC_interpolation.cc" #include "TSC_interpolation.cc" #include "SPH_interpolation.cc" @@ -270,6 +271,8 @@ void interpolate(vector *allParticles, { if ( userOptions.DTFE ) DTFE_interpolation( allParticles, samples, userOptions, uQuantities, aQuantities ); + else if ( userOptions.NGP ) + NGP_interpolation( allParticles, samples, userOptions, aQuantities ); else if ( userOptions.CIC ) CIC_interpolation( allParticles, samples, userOptions, aQuantities ); else if ( userOptions.TSC ) diff --git a/src/NGP_interpolation.cc b/src/NGP_interpolation.cc new file mode 100644 index 0000000..ac8ad54 --- /dev/null +++ b/src/NGP_interpolation.cc @@ -0,0 +1,177 @@ +/* + * Copyright (c) 2011 Marius Cautun + * + * Kapteyn Astronomical Institute + * University of Groningen, the Netherlands + * + * + * This program 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. + * + * This program 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 this program. If not, see . + * + */ + + +#include +#include +#include +#include + +#include "define.h" +#include "particle_data.h" +#include "quantities.h" +#include "box.h" +#include "user_options.h" +#include "message.h" + +using namespace std; +typedef vector::iterator vectorIterator; + + +void NGP_interpolation_regular_grid(vector &particles, + User_options &userOptions, + Quantities *q); + + + + +/* This function interpolates the density and velocity to grid using the NGP (cloud in cell) method. */ +void NGP_interpolation(vector *particles, + vector &samples, + User_options &userOptions, + Quantities *q) +{ + if ( samples.empty() and not userOptions.redshiftConeOn ) // interpolate using a regular cubic grid + NGP_interpolation_regular_grid( *particles, userOptions, q ); + else + throwError( "The NGP method can interpolate the fields only on a regular rectangular grid. No NGP interpolation methods are implemented for redshift cone coordinates or for user defined sample points." ); + particles->clear(); +} + + + + + +/* This function uses the NGP method to interpolate quantities to a grid. It interpolates only the density and the velocity. +NOTE: It does not interpolate the velocity to the grid, but in fact the momentum. The velocity is than obtained as the momentum in the cell divided by the mass in the grid cell. +*/ +void NGP_interpolation_regular_grid(vector &particles, + User_options &userOptions, + Quantities *q) +{ + size_t const *nGrid = &(userOptions.gridSize[0]); + MESSAGE::Message message( userOptions.verboseLevel ); + message << "\nInterpolating the fields to the grid using the NGP method. The interpolation takes place inside the box of coordinates " << userOptions.region.print() + << " on a " << MESSAGE::printElements( nGrid, NO_DIM, "*" ) << " grid ... " << MESSAGE::Flush; + boost::timer t; + t.restart(); + + + // allocate memory for the results + size_t const reserveSize = (NO_DIM==2) ? nGrid[0]*nGrid[1] : nGrid[0]*nGrid[1]*nGrid[2]; + q->density.assign( reserveSize, Real(0.) ); + if ( userOptions.aField.velocity ) + q->velocity.assign( reserveSize, Pvector::zero() ); + + + // get the grid spacing + Box box = userOptions.region; + Real dx[NO_DIM]; + for (int i=0; iposition(j)-box[2*j])/dx[j] ) ); + if ( cell[j]<0 or cell[j]>=nGrid[j] ) + validCell = false; + } + if ( not validCell ) continue; +#if NO_DIM==2 + int index = cell[0] * nGrid[1] + cell[1]; +#elif NO_DIM==3 + int index = cell[0] * nGrid[1]*nGrid[2] + cell[1] * nGrid[2] + cell[2]; +#endif + q->density[index] += it->weight(); + if ( userOptions.aField.velocity ) + q->velocity[index] += it->velocity() * it->weight(); + } + + + // divide the momentum by the mass in the cells + if ( userOptions.aField.velocity ) + for (size_t i=0; idensity[i]!=Real(0.) ) + q->velocity[i] /= q->density[i]; + else + q->velocity[i] = Pvector::zero(); + + // normalize the density to average background density + if ( userOptions.aField.density ) + { + Real factor = Real( q->density.size() ) / box.volume() / userOptions.averageDensity; + for (vector::iterator it=q->density.begin(); it!=q->density.end(); ++it) + (*it) *= factor; + } + else + q->density.clear(); + + message << "Done.\n" << MESSAGE::Flush; + printElapsedTime( &t, &userOptions, "NGP interpolation" ); +} + + +/* This function counts how many particles are in each cell of a grid using the NGP method. */ +void NGP_particle_count(vector &particles, + size_t const *nGrid, + Box box, + vector *counts) +{ + // allocate memory for the results + size_t const reserveSize = (NO_DIM==2) ? nGrid[0]*nGrid[1] : nGrid[0]*nGrid[1]*nGrid[2]; + counts->assign( reserveSize, int(0) ); + + // get the grid spacing + Real dx[NO_DIM]; + for (int i=0; iposition(j)-box[2*j])/dx[j] ) ); + if ( cell[j]<0 or cell[j]>=nGrid[j] ) + validCell = false; + } + if ( not validCell ) continue; +#if NO_DIM==2 + int index = cell[0] * nGrid[1] + cell[1]; +#elif NO_DIM==3 + int index = cell[0] * nGrid[1]*nGrid[2] + cell[1] * nGrid[2] + cell[2]; +#endif + (*counts)[index] += 1; + } +} + + + diff --git a/src/subpartition.h b/src/subpartition.h index 29d1a1a..e5da880 100644 --- a/src/subpartition.h +++ b/src/subpartition.h @@ -154,7 +154,7 @@ void insertParticlesInBox(std::vector &p, } -void CIC_particle_count(std::vector &particles, +void NGP_particle_count(std::vector &particles, size_t const *nGrid, Box box, std::vector *counts); @@ -245,7 +245,7 @@ void optimalPartitionSplit(std::vector &particles, // find how many particles are in each cell of the grid std::vector counts; - CIC_particle_count( particles, grid, userOptions.region, &counts); + NGP_particle_count( particles, grid, userOptions.region, &counts); // split along the x-direction std::vector xSplitIndices; xSplitIndices.assign( partition[0]+1, size_t(0) ); diff --git a/src/user_options.cc b/src/user_options.cc index 40cb132..f9d4f67 100644 --- a/src/user_options.cc +++ b/src/user_options.cc @@ -74,6 +74,7 @@ User_options::User_options() // additional options DTFE = true; + NGP = false; CIC = false; TSC = false; SPH = false; @@ -187,6 +188,7 @@ void User_options::addOptions(po::options_description &allOptions, po::options_description additionalOptions("Additional options"); additionalOptions.add_options() ("config,c", po::value(&(this->configFilename)), "supply all/part of the program options in a configuration file. The option syntax is the same as at the command line. Can insert comments using the '#' symbol (everything after this symbol until the end of the line will be considered a comment)." ) + ("NGP", "choose the NGP (nearest grid point) as the grid interpolation method instead of DTFE. This method is available only for: density and velocity fields.") ("CIC", "choose the CIC (Cloud In Cell) as the grid interpolation method instead of DTFE. This method is available only for: density and velocity fields.") ("TSC", "choose the TSC (Triangular Shape Cloud) as the grid interpolation method instead of DTFE. This method is available only for: density and velocity fields.") ("SPH", po::value(&(this->SPH_neighbors)), "choose the SPH (Smoothed Particle Hydrodynamics) as the grid interpolation method instead of DTFE. This method is available only for: density, velocity and scalar fields.") @@ -323,6 +325,7 @@ void User_options::shortHelp( char *progName ) po::options_description additionalOptions("Additional options"); additionalOptions.add_options() ("config,c", po::value< Real >(&temp), "name the configuration file from which to read the program options." ) + ("NGP", "choose NGP (Nearest Grid Point) for grid interpolation.") ("CIC", "choose CIC (Cloud In Cell) for grid interpolation.") ("TSC", "choose TSC (Triangular Shape Cloud) for grid interpolation.") ("SPH", po::value< Real >(&temp), "choose SPH for grid interpolation (argument = number nearest neighbors).") @@ -387,7 +390,8 @@ void User_options::printOptions() if ( this->uField.scalar ) uField += " scalar,"; if ( this->uField.scalar_gradient ) uField += " scalar gradient,"; if ( this->uField.triangulation ) uField += " triangulation,"; - if ( not this->uField.selected() and this->CIC ) uField = "none since CIC interpolation"; + if ( not this->uField.selected() and this->NGP ) uField = "none since NGP interpolation"; + else if ( not this->uField.selected() and this->CIC ) uField = "none since CIC interpolation"; else if ( not this->uField.selected() and this->TSC ) uField = "none since TSC interpolation"; else if ( not this->uField.selected() and this->SPH ) uField = "none since SPH interpolation"; else if ( not this->uField.selected() ) uField = "none"; @@ -406,6 +410,7 @@ void User_options::printOptions() std::string interpolationMethod = "DTFE"; if ( this->CIC ) interpolationMethod = "CIC"; + else if ( this->NGP ) interpolationMethod = "NGP"; else if ( this->TSC ) interpolationMethod = "TSC"; else if ( this->SPH ) { @@ -589,6 +594,9 @@ void User_options::readOptions(int argc, char *argv[], bool getFileNames, bool s conflicting_options(vm, "SPH", "TSC"); conflicting_options(vm, "SPH", "CIC"); conflicting_options(vm, "CIC", "TSC"); + conflicting_options(vm, "NGP", "CIC"); + conflicting_options(vm, "NGP", "TSC"); + conflicting_options(vm, "NGP", "SPH"); conflicting_options(vm, "region", "regionMpc"); conflicting_options(vm, "padding", "paddingMpc"); @@ -772,6 +780,16 @@ void User_options::readOptions(int argc, char *argv[], bool getFileNames, bool s // read the additional options + if ( vm.count("NGP") ) + { + this->NGP = true; + this->DTFE = false; + } + if ( vm.count("CIC") ) + { + this->CIC = true; + this->DTFE = false; + } if ( vm.count("TSC") ) { this->TSC = true; @@ -834,7 +852,7 @@ void User_options::readOptions(int argc, char *argv[], bool getFileNames, bool s #endif // some special settings only for the TSC or SPH methods - if ( this->CIC or this->TSC or this->SPH ) + if ( this->NGP or this->CIC or this->TSC or this->SPH ) { // for CIC, TSC and SPH only one type of fields are available - the averaged ones if ( this->uField.density ) this->aField.density = true; @@ -851,17 +869,17 @@ void User_options::readOptions(int argc, char *argv[], bool getFileNames, bool s if ( this->aField.velocity_gradient or this->aField.selectedVelocityDerivatives() ) { MESSAGE::Warning warning( verboseLevel ); - warning << "The CIC, TSC or SPH grid interpolation methods do not have implemented a method for computing the velocity gradient. The velocity gradient will not be computed!" << MESSAGE::EndWarning; + warning << "The NGP, CIC, TSC or SPH grid interpolation methods do not have implemented a method for computing the velocity gradient. The velocity gradient will not be computed!" << MESSAGE::EndWarning; this->aField.velocity_gradient = false; this->aField.deselectVelocityDerivatives(); // switches off the divergence, shear and vorticity computations } if ( this->aField.scalar_gradient ) { MESSAGE::Warning warning( verboseLevel ); - warning << "The CIC, TSC or SPH grid interpolation methods do not have implemented a method for computing the velocity gradient. The velocity gradient will not be computed!" << MESSAGE::EndWarning; + warning << "The NGP, CIC, TSC or SPH grid interpolation methods do not have implemented a method for computing the velocity gradient. The velocity gradient will not be computed!" << MESSAGE::EndWarning; this->aField.scalar_gradient = false; } - if ( (this->CIC or this->TSC) and this->aField.scalar ) + if ( (this->NGP or this->CIC or this->TSC) and this->aField.scalar ) { MESSAGE::Warning warning( verboseLevel ); warning << "The CIC and TSC grid interpolation method does not have implemented a method for computing the scalar fields on the grid. The scalar field cannot be computed!" << MESSAGE::EndWarning; @@ -1007,14 +1025,14 @@ void User_options::updateEntries(size_t const noTotalParticles, // switch gradient computations off for the TSC and SPH methods - if ( CIC or TSC or SPH ) + if ( NGP or CIC or TSC or SPH ) { if ( aField.velocity_gradient or aField.selectedVelocityDerivatives() ) - throwError( "The CIC, TSC and SPH grid interpolation methods do not have implemented a method for computing the velocity gradient. The velocity gradient cannot be computed!"); + throwError( "The NGP, CIC, TSC and SPH grid interpolation methods do not have implemented a method for computing the velocity gradient. The velocity gradient cannot be computed!"); if ( aField.scalar_gradient ) - throwError( "The CIC, TSC and SPH grid interpolation methods do not have implemented a method for computing the scalar fields gradients. The scalar gradient cannot be computed!"); + throwError( "The NGP, CIC, TSC and SPH grid interpolation methods do not have implemented a method for computing the scalar fields gradients. The scalar gradient cannot be computed!"); if ( (CIC or TSC) and aField.scalar ) - throwError( "The CIC and TSC grid interpolation method does not have implemented a method for computing the scalar fields on the grid. The scalar field cannot be computed!"); + throwError( "The NGP, CIC and TSC grid interpolation method does not have implemented a method for computing the scalar fields on the grid. The scalar field cannot be computed!"); } @@ -1047,7 +1065,7 @@ void User_options::updatePadding(size_t const noTotalParticles) if ( this->DTFE or this->SPH ) //padding for the DTFE and SPH methods for (int i=0; i<2*NO_DIM; ++i) paddingLength[i] = paddingParticles / Real( pow(noTotalParticles,1./NO_DIM) ) * fullBoxLength[i%2]; - else if ( this->TSC ) + else if ( this->TSC or this->NGP ) for (int i=0; i<2*NO_DIM; ++i) paddingLength[i] = fullBoxLength[i%2] / gridSize[i%2]; else if ( this->CIC ) diff --git a/src/user_options.h b/src/user_options.h index 04738b3..6275c27 100644 --- a/src/user_options.h +++ b/src/user_options.h @@ -152,6 +152,7 @@ struct User_options // additional options std::string configFilename;// the program options can be supplied in a file too. This keeps track of the file name where the options where supplied. + bool NGP; // true if to use CIC grid interpolation instead of DTFE bool CIC; // true if to use CIC grid interpolation instead of DTFE bool TSC; // true if to use TSC grid interpolation instead of DTFE bool SPH; // true if to use SPH grid interpolation instead of DTFE