From 08d8298428c9839f5dcaf050e6eb2950674179a0 Mon Sep 17 00:00:00 2001
From: Gareth Tribello <gareth.tribello@gmail.com>
Date: Sat, 30 Jan 2016 12:22:17 +0000
Subject: [PATCH] Added first go at faster wcsurface stuff
Now there is an option to update grids only in the region around where
the contour was found in FindContour
---
regtest/multicolvar/rt-wcsurface/plumed.dat | 9 ------
src/gridtools/ActionWithInputGrid.cpp | 6 ++++
src/gridtools/ActionWithInputGrid.h | 1 +
src/gridtools/FindContour.cpp | 33 +++++++++++++++++++--
src/gridtools/GridVessel.cpp | 9 ++++--
src/gridtools/GridVessel.h | 22 ++++++++++----
src/gridtools/HistogramOnGrid.cpp | 1 +
7 files changed, 63 insertions(+), 18 deletions(-)
delete mode 100644 regtest/multicolvar/rt-wcsurface/plumed.dat
diff --git a/regtest/multicolvar/rt-wcsurface/plumed.dat b/regtest/multicolvar/rt-wcsurface/plumed.dat
deleted file mode 100644
index 452ac0e42..000000000
--- a/regtest/multicolvar/rt-wcsurface/plumed.dat
+++ /dev/null
@@ -1,9 +0,0 @@
-UNITS NATURAL
-
-fcc: FCCUBIC SPECIES=1-5184 SWITCH={CUBIC D_0=1.2 D_MAX=1.5} ALPHA=27
-
-DUMPMULTICOLVAR DATA=fcc ORIGIN=1 FILE=atoms.xyz PRECISION=4
-dens: MULTICOLVARDENS DATA=fcc ORIGIN=1 DIR=xyz NBINS=14,14,28 BANDWIDTH=1.0,1.0,1.0
-
-PRINT_CUBE GRID=dens STRIDE=1 FILE=dens.cube
-FIND_CONTOUR GRID=dens STRIDE=1 CONTOUR=0.42 FILE=surface.xyz PRECISION=4
diff --git a/src/gridtools/ActionWithInputGrid.cpp b/src/gridtools/ActionWithInputGrid.cpp
index 7596d020b..88bc3ee69 100644
--- a/src/gridtools/ActionWithInputGrid.cpp
+++ b/src/gridtools/ActionWithInputGrid.cpp
@@ -78,6 +78,12 @@ void ActionWithInputGrid::setAnalysisStride( const bool& use_all, const unsigned
void ActionWithInputGrid::update(){
if( unormalised ) norm = 1.0;
else norm=1.0/mygrid->getNorm();
+
+ if( checkAllActive() ){
+ for(unsigned i=0;i<mygrid->getNumberOfPoints();++i){
+ if( mygrid->inactive(i) ) error("if FIND_CONTOUR is used with BUFFER option then other actions cannot be performed with grid");
+ }
+ }
performOperationsWithGrid( true );
}
diff --git a/src/gridtools/ActionWithInputGrid.h b/src/gridtools/ActionWithInputGrid.h
index 9b5251641..0246be6ae 100644
--- a/src/gridtools/ActionWithInputGrid.h
+++ b/src/gridtools/ActionWithInputGrid.h
@@ -47,6 +47,7 @@ public:
void apply(){}
void update();
void runFinalJobs();
+ virtual bool checkAllActive() const { return true; }
virtual void performOperationsWithGrid( const bool& from_update )=0;
void setAnalysisStride( const bool& use_all, const unsigned& astride );
};
diff --git a/src/gridtools/FindContour.cpp b/src/gridtools/FindContour.cpp
index 0a1c135ed..1ba4c7bce 100644
--- a/src/gridtools/FindContour.cpp
+++ b/src/gridtools/FindContour.cpp
@@ -33,6 +33,7 @@ class FindContour : public ActionWithInputGrid {
private:
OFile of;
double lenunit;
+ unsigned gbuffer;
std::string fmt_xyz;
unsigned mycomp;
double contour;
@@ -40,6 +41,7 @@ private:
public:
static void registerKeywords( Keywords& keys );
explicit FindContour(const ActionOptions&ao);
+ bool checkAllActive() const { return gbuffer==0; }
void performOperationsWithGrid( const bool& from_update );
double getDifferenceFromContour( const std::vector<double>& x, std::vector<double>& der );
unsigned getNumberOfDerivatives(){ return 0; }
@@ -61,6 +63,7 @@ void FindContour::registerKeywords( Keywords& keys ){
"is thre dimensional specifying SEARCHDIR=1,2 means that only points found along grid lines parallel to the "
"first and second axis of the grid will be looked for. The code will not search for the dividing surface "
"along the third axis of the grid");
+ keys.add("compulsory","BUFFER","0","number of buffer grid points around location where grid was found on last step. If this is zero the full grid is calculated on each step");
keys.add("optional","COMPONENT","if your input is a vector field use this to specifiy the component of the input vector field for which you wish to find the contour");
keys.add("optional", "PRECISION","The number of digits in trajectory file");
}
@@ -79,8 +82,9 @@ nosearch_dirs( mygrid->getDimension() )
}
if( mygrid->noDerivatives() ) error("cannot find contours if input grid has no derivatives");
- parse("CONTOUR",contour);
+ parse("CONTOUR",contour); parse("BUFFER",gbuffer);
log.printf(" calculating dividing surface along which function equals %f \n", contour);
+ if( gbuffer>0 ) log.printf(" after first step a subset of only %d grid points around where the countour was found will be checked\n",gbuffer);
std::string searchdir_str; parse("SEARCHDIR",searchdir_str);
if( searchdir_str=="all" ){
@@ -149,10 +153,28 @@ void FindContour::performOperationsWithGrid( const bool& from_update ){
}
// Run over whole grid
+ std::vector<bool> active( mygrid->getNumberOfPoints(), false );
+ std::vector<unsigned> neighbours; unsigned num_neighbours;
+ std::vector<unsigned> gbuffer_vec( mygrid->getDimension(), gbuffer );
+ std::vector<unsigned> ugrid_indices( mygrid->getDimension() );
unsigned npoints=0; RootFindingBase<FindContour> mymin( this );
for(unsigned i=0;i<mygrid->getNumberOfPoints();++i){
// Get the index of the current grid point
mygrid->getIndices( i, ind );
+ // Ensure inactive grid points are ignored
+ bool cycle=false;
+ if( mygrid->inactive(i) ) continue ;
+ for(unsigned j=0;j<mygrid->getDimension();++j) ugrid_indices[j]=ind[j];
+ for(unsigned j=0;j<mygrid->getDimension();++j){
+ if( mygrid->isPeriodic(j) ) ugrid_indices[j]=(ugrid_indices[j]+1)%nbin[j];
+ else if( (ugrid_indices[j]+1)==nbin[j] ) continue ;
+ else ugrid_indices[j]+=1;
+ // Now check the grid is active
+ if( mygrid->inactive( mygrid->getIndex( ugrid_indices ) ) ){ cycle=true; break; }
+ // And reset ugrid_indices
+ ugrid_indices[j]=ind[j];
+ }
+ if( cycle ) continue ;
// Get the value of a point on the grid
double val1=getGridElement( i, mycomp*(mygrid->getDimension()+1) ) - contour;
@@ -173,13 +195,20 @@ void FindContour::performOperationsWithGrid( const bool& from_update ){
direction[j]=0.999999999*dx[j];
// And do proper search for contour point
mymin.linesearch( direction, contour_points[npoints], &FindContour::getDifferenceFromContour );
- direction[j]=0.0; npoints++;
+ direction[j]=0.0;
+ // For next run through activate buffer region around this grid point
+ if( gbuffer>0 ){
+ mygrid->getNeighbors( contour_points[npoints], gbuffer_vec, num_neighbours, neighbours );
+ for(unsigned n=0;n<num_neighbours;++n) active[ neighbours[n] ]=true;
+ }
+ npoints++;
}
if( mygrid->isPeriodic(j) && edge ){ edge=false; ind[j]=nbin[j]-1; }
else ind[j]-=1;
}
}
+ if( gbuffer>0 ) mygrid->activateThesePoints( active );
// Clear the grid ready for next time
if( from_update ) mygrid->reset();
diff --git a/src/gridtools/GridVessel.cpp b/src/gridtools/GridVessel.cpp
index 6c3b498b0..36b3e69a2 100644
--- a/src/gridtools/GridVessel.cpp
+++ b/src/gridtools/GridVessel.cpp
@@ -120,7 +120,7 @@ std::string GridVessel::description(){
void GridVessel::resize(){
plumed_massert( nper>0, "Number of datapoints at each grid point has not been set");
- resizeBuffer( npoints*nper ); data.resize( npoints*nper, 0 );
+ resizeBuffer( npoints*nper ); data.resize( npoints*nper, 0 ); active.resize( npoints, true );
}
unsigned GridVessel::getIndex( const std::vector<unsigned>& indices ) const {
@@ -189,7 +189,7 @@ void GridVessel::getSplineNeighbors( const unsigned& mybox, std::vector<unsigned
}
double GridVessel::getGridElement( const unsigned& ipoint, const unsigned& jelement ) const {
- plumed_dbg_assert( bounds_set && ipoint<npoints && jelement<nper );
+ plumed_assert( bounds_set && ipoint<npoints && jelement<nper && active[ipoint] );
return data[ nper*ipoint + jelement ];
}
@@ -344,6 +344,11 @@ double GridVessel::getValueAndDerivatives( const std::vector<double>& x, const u
return value;
}
+void GridVessel::activateThesePoints( const std::vector<bool>& to_activate ){
+ plumed_dbg_assert( to_activate.size()==npoints );
+ for(unsigned i=0;i<npoints;++i) active[i]=to_activate[i];
+}
+
}
}
diff --git a/src/gridtools/GridVessel.h b/src/gridtools/GridVessel.h
index cdc74eac6..2e17bcbe4 100644
--- a/src/gridtools/GridVessel.h
+++ b/src/gridtools/GridVessel.h
@@ -50,8 +50,6 @@ private:
std::vector<unsigned> stride;
/// The number of bins in each grid direction
std::vector<unsigned> nbin;
-/// Flatten the grid and get the grid index for a point
- unsigned getIndex( const std::vector<unsigned>& indices ) const ;
/// The grid point that was requested last by getGridPointCoordinates
unsigned currentGridPoint;
protected:
@@ -68,11 +66,10 @@ protected:
std::vector<double> dx;
/// The dimensionality of the grid
unsigned dimension;
+/// Which grid points are we actively accumulating
+ std::vector<bool> active;
/// The grid with all the data values on it
std::vector<double> data;
-/// Get the set of points neighouring a particular location in space
- void getNeighbors( const std::vector<double>& pp, const std::vector<unsigned>& nneigh,
- unsigned& num_neighbours, std::vector<unsigned>& neighbors ) const ;
/// Get the indices of a particular point
void getIndices( const std::vector<double>& point, std::vector<unsigned>& indices ) const ;
/// Convert a point in space the the correspoinding grid point
@@ -90,6 +87,8 @@ public:
std::string description();
/// Convert an index into indices
void convertIndexToIndices( const unsigned& index, const std::vector<unsigned>& nnbin, std::vector<unsigned>& indices ) const ;
+/// Flatten the grid and get the grid index for a point
+ unsigned getIndex( const std::vector<unsigned>& indices ) const ;
/// Get the indices fof a point
void getIndices( const unsigned& index, std::vector<unsigned>& indices ) const ;
@@ -127,6 +126,9 @@ public:
double getCellVolume() const ;
/// Get the value of the ith grid element
double getGridElement( const unsigned&, const unsigned& ) const ;
+/// Get the set of points neighouring a particular location in space
+ void getNeighbors( const std::vector<double>& pp, const std::vector<unsigned>& nneigh,
+ unsigned& num_neighbours, std::vector<unsigned>& neighbors ) const ;
/// Get the points neighboring a particular spline point
void getSplineNeighbors( const unsigned& mybox, std::vector<unsigned>& mysneigh ) const ;
/// Get the spacing between grid points
@@ -152,6 +154,10 @@ public:
double getValueAndDerivatives( const std::vector<double>& x, const unsigned& ind, std::vector<double>& der ) const ;
/// Was the grid cleared on the last step
bool wasreset() const ;
+/// Deactivate all the grid points
+ void activateThesePoints( const std::vector<bool>& to_activate );
+/// Is this point active
+ bool inactive( const unsigned& ip ) const ;
};
inline
@@ -216,6 +222,12 @@ bool GridVessel::wasreset() const {
return wascleared;
}
+inline
+bool GridVessel::inactive( const unsigned& ip ) const {
+ plumed_dbg_assert( ip<npoints );
+ return !active[ip];
+}
+
}
}
#endif
diff --git a/src/gridtools/HistogramOnGrid.cpp b/src/gridtools/HistogramOnGrid.cpp
index edb89d003..a973ef906 100644
--- a/src/gridtools/HistogramOnGrid.cpp
+++ b/src/gridtools/HistogramOnGrid.cpp
@@ -85,6 +85,7 @@ bool HistogramOnGrid::calculate( const unsigned& current, MultiValue& myvals, st
double newval; std::vector<double> der( dimension );
for(unsigned i=0;i<num_neigh;++i){
unsigned ineigh=neighbors[i];
+ if( inactive( ineigh ) ) continue ;
getGridPointCoordinates( ineigh, xx );
for(unsigned j=0;j<dimension;++j) vv[j]->set(xx[j]);
newval = kernel.evaluate( vv, der, true );
--
GitLab