LCOV - code coverage report
Current view: top level - vesselbase - ActionWithVessel.cpp (source / functions) Hit Total Coverage
Test: plumed test coverage Lines: 196 208 94.2 %
Date: 2020-11-18 11:20:57 Functions: 23 29 79.3 %

          Line data    Source code
       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2012-2019 The plumed team
       3             :    (see the PEOPLE file at the root of the distribution for a list of names)
       4             : 
       5             :    See http://www.plumed.org for more information.
       6             : 
       7             :    This file is part of plumed, version 2.
       8             : 
       9             :    plumed is free software: you can redistribute it and/or modify
      10             :    it under the terms of the GNU Lesser General Public License as published by
      11             :    the Free Software Foundation, either version 3 of the License, or
      12             :    (at your option) any later version.
      13             : 
      14             :    plumed is distributed in the hope that it will be useful,
      15             :    but WITHOUT ANY WARRANTY; without even the implied warranty of
      16             :    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      17             :    GNU Lesser General Public License for more details.
      18             : 
      19             :    You should have received a copy of the GNU Lesser General Public License
      20             :    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : #include "ActionWithVessel.h"
      23             : #include "tools/Communicator.h"
      24             : #include "Vessel.h"
      25             : #include "ShortcutVessel.h"
      26             : #include "StoreDataVessel.h"
      27             : #include "VesselRegister.h"
      28             : #include "BridgeVessel.h"
      29             : #include "FunctionVessel.h"
      30             : #include "StoreDataVessel.h"
      31             : #include "tools/OpenMP.h"
      32             : #include "tools/Stopwatch.h"
      33             : 
      34             : using namespace std;
      35             : namespace PLMD {
      36             : namespace vesselbase {
      37             : 
      38         513 : void ActionWithVessel::registerKeywords(Keywords& keys) {
      39        2052 :   keys.add("hidden","TOL","this keyword can be used to speed up your calculation. When accumulating sums in which the individual "
      40             :            "terms are numbers inbetween zero and one it is assumed that terms less than a certain tolerance "
      41             :            "make only a small contribution to the sum.  They can thus be safely ignored as can the the derivatives "
      42             :            "wrt these small quantities.");
      43        2052 :   keys.add("hidden","MAXDERIVATIVES","The maximum number of derivatives that can be used when storing data.  This controls when "
      44             :            "we have to start using lowmem");
      45        1539 :   keys.addFlag("SERIAL",false,"do the calculation in serial.  Do not parallelize");
      46        1539 :   keys.addFlag("LOWMEM",false,"lower the memory requirements");
      47        1539 :   keys.addFlag("TIMINGS",false,"output information on the timings of the various parts of the calculation");
      48        1539 :   keys.reserveFlag("HIGHMEM",false,"use a more memory intensive version of this collective variable");
      49         513 :   keys.add( vesselRegister().getKeywords() );
      50         513 : }
      51             : 
      52         422 : ActionWithVessel::ActionWithVessel(const ActionOptions&ao):
      53             :   Action(ao),
      54             :   serial(false),
      55             :   lowmem(false),
      56             :   noderiv(true),
      57             :   actionIsBridged(false),
      58             :   nactive_tasks(0),
      59         422 :   stopwatch(*new Stopwatch),
      60             :   dertime_can_be_off(false),
      61             :   dertime(true),
      62             :   contributorsAreUnlocked(false),
      63             :   weightHasDerivatives(false),
      64        1266 :   mydata(NULL)
      65             : {
      66         844 :   maxderivatives=309; parse("MAXDERIVATIVES",maxderivatives);
      67        1264 :   if( keywords.exists("SERIAL") ) parseFlag("SERIAL",serial);
      68           2 :   else serial=true;
      69         422 :   if(serial)log.printf("  doing calculation in serial\n");
      70         844 :   if( keywords.exists("LOWMEM") ) {
      71         786 :     plumed_assert( !keywords.exists("HIGHMEM") );
      72         786 :     parseFlag("LOWMEM",lowmem);
      73         393 :     if(lowmem) {
      74          26 :       log.printf("  lowering memory requirements\n");
      75          26 :       dertime_can_be_off=true;
      76             :     }
      77             :   }
      78         844 :   if( keywords.exists("HIGHMEM") ) {
      79          42 :     plumed_assert( !keywords.exists("LOWMEM") );
      80          42 :     bool highmem; parseFlag("HIGHMEM",highmem);
      81          21 :     lowmem=!highmem;
      82          21 :     if(!lowmem) log.printf("  increasing the memory requirements\n");
      83             :   }
      84         422 :   tolerance=nl_tolerance=epsilon;
      85        1261 :   if( keywords.exists("TOL") ) parse("TOL",tolerance);
      86         422 :   if( tolerance>epsilon) {
      87           6 :     log.printf(" Ignoring contributions less than %f \n",tolerance);
      88             :   }
      89         844 :   parseFlag("TIMINGS",timers);
      90         422 :   stopwatch.start(); stopwatch.pause();
      91         422 : }
      92             : 
      93         422 : ActionWithVessel::~ActionWithVessel() {
      94        2416 :   for(unsigned i=0; i<functions.size(); ++i) delete functions[i];
      95         422 :   stopwatch.start(); stopwatch.stop();
      96         422 :   if(timers) {
      97           0 :     log.printf("timings for action %s with label %s \n", getName().c_str(), getLabel().c_str() );
      98           0 :     log<<stopwatch;
      99             :   }
     100         844 :   delete &stopwatch;
     101         422 : }
     102             : 
     103         328 : void ActionWithVessel::addVessel( const std::string& name, const std::string& input, const int numlab ) {
     104         984 :   VesselOptions da(name,"",numlab,input,this);
     105         984 :   Vessel* vv=vesselRegister().create(name,da);
     106         328 :   FunctionVessel* fv=dynamic_cast<FunctionVessel*>(vv);
     107         328 :   if( fv ) {
     108         303 :     std::string mylabel=Vessel::transformName( name );
     109         303 :     plumed_massert( keywords.outputComponentExists(mylabel,false), "a description of the value calculated by vessel " + name + " has not been added to the manual");
     110             :   }
     111         328 :   addVessel(vv);
     112         328 : }
     113             : 
     114         488 : void ActionWithVessel::addVessel( Vessel* vv ) {
     115         488 :   ShortcutVessel* sv=dynamic_cast<ShortcutVessel*>(vv);
     116         488 :   if(!sv) { vv->checkRead(); functions.push_back(vv); }
     117          10 :   else { delete sv; return; }
     118             : 
     119         956 :   StoreDataVessel* mm=dynamic_cast<StoreDataVessel*>( vv );
     120         478 :   if( mydata && mm ) error("cannot have more than one StoreDataVessel in one action");
     121         478 :   else if( mm ) mydata=mm;
     122         359 :   else dertime_can_be_off=false;
     123             : }
     124             : 
     125          46 : BridgeVessel* ActionWithVessel::addBridgingVessel( ActionWithVessel* tome ) {
     126         230 :   VesselOptions da("","",0,"",this);
     127          46 :   BridgeVessel* bv=new BridgeVessel(da);
     128          46 :   bv->setOutputAction( tome );
     129          46 :   tome->actionIsBridged=true; dertime_can_be_off=false;
     130          92 :   functions.push_back( dynamic_cast<Vessel*>(bv) );
     131          46 :   resizeFunctions();
     132          46 :   return bv;
     133             : }
     134             : 
     135         171 : StoreDataVessel* ActionWithVessel::buildDataStashes( ActionWithVessel* actionThatUses ) {
     136         171 :   if(mydata) {
     137          81 :     if( actionThatUses ) mydata->addActionThatUses( actionThatUses );
     138          81 :     return mydata;
     139             :   }
     140             : 
     141         450 :   VesselOptions da("","",0,"",this);
     142          90 :   StoreDataVessel* mm=new StoreDataVessel(da);
     143          90 :   if( actionThatUses ) mm->addActionThatUses( actionThatUses );
     144          90 :   addVessel(mm);
     145             : 
     146             :   // Make sure resizing of vessels is done
     147          90 :   resizeFunctions();
     148             : 
     149          90 :   return mydata;
     150             : }
     151             : 
     152    12650372 : void ActionWithVessel::addTaskToList( const unsigned& taskCode ) {
     153    25300744 :   fullTaskList.push_back( taskCode ); taskFlags.push_back(0);
     154    12650372 :   plumed_assert( fullTaskList.size()==taskFlags.size() );
     155    12650372 : }
     156             : 
     157         391 : void ActionWithVessel::readVesselKeywords() {
     158             :   // Set maxderivatives if it is too big
     159         391 :   if( maxderivatives>getNumberOfDerivatives() ) maxderivatives=getNumberOfDerivatives();
     160             : 
     161             :   // Loop over all keywords find the vessels and create appropriate functions
     162       16849 :   for(unsigned i=0; i<keywords.size(); ++i) {
     163       16458 :     std::string thiskey,input; thiskey=keywords.getKeyword(i);
     164             :     // Check if this is a key for a vessel
     165       24687 :     if( vesselRegister().check(thiskey) ) {
     166        5002 :       plumed_assert( keywords.style(thiskey,"vessel") );
     167        2501 :       bool dothis=false; parseFlag(thiskey,dothis);
     168        2501 :       if(dothis) addVessel( thiskey, input );
     169             : 
     170        2501 :       parse(thiskey,input);
     171        2501 :       if(input.size()!=0) {
     172         115 :         addVessel( thiskey, input );
     173             :       } else {
     174          22 :         for(unsigned i=1;; ++i) {
     175        2408 :           if( !parseNumbered(thiskey,i,input) ) break;
     176          22 :           std::string ss; Tools::convert(i,ss);
     177          22 :           addVessel( thiskey, input, i );
     178             :           input.clear();
     179          22 :         }
     180             :       }
     181             :     }
     182             :   }
     183             : 
     184             :   // Make sure all vessels have had been resized at start
     185         391 :   if( functions.size()>0 ) resizeFunctions();
     186         391 : }
     187             : 
     188        1274 : void ActionWithVessel::resizeFunctions() {
     189        8518 :   for(unsigned i=0; i<functions.size(); ++i) functions[i]->resize();
     190        1274 : }
     191             : 
     192         765 : void ActionWithVessel::needsDerivatives() {
     193             :   // Turn on the derivatives and resize
     194         765 :   noderiv=false; resizeFunctions();
     195             :   // Setting contributors unlocked here ensures that link cells are ignored
     196         765 :   contributorsAreUnlocked=true; contributorsAreUnlocked=false;
     197             :   // And turn on the derivatives in all actions on which we are dependent
     198        2373 :   for(unsigned i=0; i<getDependencies().size(); ++i) {
     199         281 :     ActionWithVessel* vv=dynamic_cast<ActionWithVessel*>( getDependencies()[i] );
     200         281 :     if(vv) vv->needsDerivatives();
     201             :   }
     202         765 : }
     203             : 
     204        3542 : void ActionWithVessel::lockContributors() {
     205        3542 :   nactive_tasks = 0;
     206    55576579 :   for(unsigned i=0; i<fullTaskList.size(); ++i) {
     207    18523165 :     if( taskFlags[i]>0 ) nactive_tasks++;
     208             :   }
     209             : 
     210             :   unsigned n=0;
     211        3542 :   partialTaskList.resize( nactive_tasks );
     212        3542 :   indexOfTaskInFullList.resize( nactive_tasks );
     213    55576579 :   for(unsigned i=0; i<fullTaskList.size(); ++i) {
     214             :     // Deactivate sets inactive tasks to number not equal to zero
     215    18523165 :     if( taskFlags[i]>0 ) {
     216    11032342 :       partialTaskList[n] = fullTaskList[i];
     217     5516171 :       indexOfTaskInFullList[n]=i;
     218     5516171 :       n++;
     219             :     }
     220             :   }
     221             :   plumed_dbg_assert( n==nactive_tasks );
     222       20980 :   for(unsigned i=0; i<functions.size(); ++i) {
     223        4632 :     BridgeVessel* bb = dynamic_cast<BridgeVessel*>( functions[i] );
     224        4632 :     if( bb ) bb->copyTaskFlags();
     225             :   }
     226             :   // Resize mydata to accomodate all active tasks
     227        3542 :   if( mydata ) mydata->resize();
     228        3542 :   contributorsAreUnlocked=false;
     229        3542 : }
     230             : 
     231        3542 : void ActionWithVessel::deactivateAllTasks() {
     232        3542 :   contributorsAreUnlocked=true; nactive_tasks = 0;
     233        7084 :   taskFlags.assign(taskFlags.size(),0);
     234        3542 : }
     235             : 
     236      205651 : bool ActionWithVessel::taskIsCurrentlyActive( const unsigned& index ) const {
     237      411302 :   plumed_dbg_assert( index<taskFlags.size() ); return (taskFlags[index]>0);
     238             : }
     239             : 
     240       22917 : void ActionWithVessel::doJobsRequiredBeforeTaskList() {
     241             :   // Do any preparatory stuff for functions
     242      148833 :   for(unsigned j=0; j<functions.size(); ++j) functions[j]->prepare();
     243       22917 : }
     244             : 
     245       23430 : unsigned ActionWithVessel::getSizeOfBuffer( unsigned& bufsize ) {
     246      151887 :   for(unsigned i=0; i<functions.size(); ++i) functions[i]->setBufferStart( bufsize );
     247       23430 :   if( buffer.size()!=bufsize ) buffer.resize( bufsize );
     248       23430 :   if( mydata ) {
     249             :     unsigned dsize=mydata->getSizeOfDerivativeList();
     250        2315 :     if( der_list.size()!=dsize ) der_list.resize( dsize );
     251             :   }
     252       23430 :   return bufsize;
     253             : }
     254             : 
     255       19892 : void ActionWithVessel::runAllTasks() {
     256       39784 :   plumed_massert( !contributorsAreUnlocked && functions.size()>0, "you must have a call to readVesselKeywords somewhere" );
     257       19892 :   unsigned stride=comm.Get_size();
     258       19892 :   unsigned rank=comm.Get_rank();
     259       19892 :   if(serial) { stride=1; rank=0; }
     260             : 
     261             :   // Make sure jobs are done
     262       19892 :   if(timers) stopwatch.start("1 Prepare Tasks");
     263       19892 :   doJobsRequiredBeforeTaskList();
     264       19892 :   if(timers) stopwatch.stop("1 Prepare Tasks");
     265             : 
     266             :   // Get number of threads for OpenMP
     267       19892 :   unsigned nt=OpenMP::getNumThreads();
     268       19892 :   if( nt*stride*2>nactive_tasks || !threadSafe()) nt=1;
     269             : 
     270             :   // Get size for buffer
     271       19892 :   unsigned bsize=0, bufsize=getSizeOfBuffer( bsize );
     272             :   // Clear buffer
     273       39784 :   buffer.assign( buffer.size(), 0.0 );
     274             :   // Switch off calculation of derivatives in main loop
     275       19892 :   if( dertime_can_be_off ) dertime=false;
     276             : 
     277       19892 :   if(timers) stopwatch.start("2 Loop over tasks");
     278       49316 :   #pragma omp parallel num_threads(nt)
     279             :   {
     280             :     std::vector<double> omp_buffer;
     281       48488 :     if( nt>1 ) omp_buffer.resize( bufsize, 0.0 );
     282       58848 :     MultiValue myvals( getNumberOfQuantities(), getNumberOfDerivatives() );
     283       58848 :     MultiValue bvals( getNumberOfQuantities(), getNumberOfDerivatives() );
     284       29424 :     myvals.clearAll(); bvals.clearAll();
     285             : 
     286             :     #pragma omp for nowait schedule(dynamic)
     287       29424 :     for(unsigned i=rank; i<nactive_tasks; i+=stride) {
     288             :       // Calculate the stuff in the loop for this action
     289     1982494 :       performTask( indexOfTaskInFullList[i], partialTaskList[i], myvals );
     290             : 
     291             :       // Check for conditions that allow us to just to skip the calculation
     292             :       // the condition is that the weight of the contribution is low
     293             :       // N.B. Here weights are assumed to be between zero and one
     294     1569085 :       if( myvals.get(0)<tolerance ) {
     295             :         // Clear the derivatives
     296      577838 :         myvals.clearAll();
     297      577838 :         continue;
     298             :       }
     299             : 
     300             :       // Now calculate all the functions
     301             :       // If the contribution of this quantity is very small at neighbour list time ignore it
     302             :       // untill next neighbour list time
     303      413409 :       if( nt>1 ) {
     304      739826 :         calculateAllVessels( indexOfTaskInFullList[i], myvals, bvals, omp_buffer, der_list );
     305             :       } else {
     306       86992 :         calculateAllVessels( indexOfTaskInFullList[i], myvals, bvals, buffer, der_list );
     307             :       }
     308             : 
     309             :       // Clear the value
     310      413409 :       myvals.clearAll();
     311             :     }
     312       58848 :     #pragma omp critical
     313    38538680 :     if(nt>1) for(unsigned i=0; i<bufsize; ++i) buffer[i]+=omp_buffer[i];
     314             :   }
     315       19892 :   if(timers) stopwatch.stop("2 Loop over tasks");
     316             :   // Turn back on derivative calculation
     317       19892 :   dertime=true;
     318             : 
     319       19892 :   if(timers) stopwatch.start("3 MPI gather");
     320             :   // MPI Gather everything
     321       39784 :   if( !serial && buffer.size()>0 ) comm.Sum( buffer );
     322             :   // MPI Gather index stores
     323       19892 :   if( mydata && !lowmem && !noderiv ) {
     324         671 :     comm.Sum( der_list ); mydata->setActiveValsAndDerivatives( der_list );
     325             :   }
     326             :   // Update the elements that are makign contributions to the sum here
     327             :   // this causes problems if we do it in prepare
     328       19892 :   if(timers) stopwatch.stop("3 MPI gather");
     329             : 
     330       19892 :   if(timers) stopwatch.start("4 Finishing computations");
     331       19892 :   finishComputations( buffer );
     332       19892 :   if(timers) stopwatch.stop("4 Finishing computations");
     333       19892 : }
     334             : 
     335           0 : void ActionWithVessel::transformBridgedDerivatives( const unsigned& current, MultiValue& invals, MultiValue& outvals ) const {
     336           0 :   plumed_error();
     337             : }
     338             : 
     339      475917 : void ActionWithVessel::calculateAllVessels( const unsigned& taskCode, MultiValue& myvals, MultiValue& bvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) {
     340     2776587 :   for(unsigned j=0; j<functions.size(); ++j) {
     341             :     // Calculate returns a bool that tells us if this particular
     342             :     // quantity is contributing more than the tolerance
     343      608251 :     functions[j]->calculate( taskCode, functions[j]->transformDerivatives(taskCode, myvals, bvals), buffer, der_list );
     344      608251 :     if( !actionIsBridged ) bvals.clearAll();
     345             :   }
     346      475917 :   return;
     347             : }
     348             : 
     349       22917 : void ActionWithVessel::finishComputations( const std::vector<double>& buffer ) {
     350             :   // Set the final value of the function
     351      148833 :   for(unsigned j=0; j<functions.size(); ++j) functions[j]->finish( buffer );
     352       22917 : }
     353             : 
     354        6971 : bool ActionWithVessel::getForcesFromVessels( std::vector<double>& forcesToApply ) {
     355             : #ifndef NDEBUG
     356             :   if( forcesToApply.size()>0 ) plumed_dbg_assert( forcesToApply.size()==getNumberOfDerivatives() );
     357             : #endif
     358        6971 :   if(tmpforces.size()!=forcesToApply.size() ) tmpforces.resize( forcesToApply.size() );
     359             : 
     360       13942 :   forcesToApply.assign( forcesToApply.size(),0.0 );
     361             :   bool wasforced=false;
     362       33301 :   for(unsigned i=0; i<getNumberOfVessels(); ++i) {
     363       26330 :     if( (functions[i]->applyForce( tmpforces )) ) {
     364             :       wasforced=true;
     365     2096364 :       for(unsigned j=0; j<forcesToApply.size(); ++j) forcesToApply[j]+=tmpforces[j];
     366             :     }
     367             :   }
     368        6971 :   return wasforced;
     369             : }
     370             : 
     371           0 : void ActionWithVessel::retrieveDomain( std::string& min, std::string& max ) {
     372           0 :   plumed_merror("If your function is periodic you need to add a retrieveDomain function so that ActionWithVessel can retrieve the domain");
     373             : }
     374             : 
     375           0 : Vessel* ActionWithVessel::getVesselWithName( const std::string& mynam ) {
     376             :   int target=-1;
     377           0 :   for(unsigned i=0; i<functions.size(); ++i) {
     378           0 :     if( functions[i]->getName().find(mynam)!=std::string::npos ) {
     379           0 :       if( target<0 ) target=i;
     380           0 :       else error("found more than one " + mynam + " object in action");
     381             :     }
     382             :   }
     383           0 :   return functions[target];
     384             : }
     385             : 
     386             : }
     387        4839 : }

Generated by: LCOV version 1.13