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 509 : void ActionWithVessel::registerKeywords(Keywords& keys) {
39 : 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 509 : "wrt these small quantities.");
43 : keys.add("hidden","MAXDERIVATIVES","The maximum number of derivatives that can be used when storing data. This controls when "
44 509 : "we have to start using lowmem");
45 509 : keys.addFlag("SERIAL",false,"do the calculation in serial. Do not parallelize");
46 509 : keys.addFlag("LOWMEM",false,"lower the memory requirements");
47 509 : keys.addFlag("TIMINGS",false,"output information on the timings of the various parts of the calculation");
48 509 : keys.reserveFlag("HIGHMEM",false,"use a more memory intensive version of this collective variable");
49 509 : keys.add( vesselRegister().getKeywords() );
50 509 : }
51 :
52 418 : 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 418 : stopwatch(*new Stopwatch),
60 : dertime_can_be_off(false),
61 : dertime(true),
62 : contributorsAreUnlocked(false),
63 : weightHasDerivatives(false),
64 836 : mydata(NULL)
65 : {
66 418 : maxderivatives=309; parse("MAXDERIVATIVES",maxderivatives);
67 418 : if( keywords.exists("SERIAL") ) parseFlag("SERIAL",serial);
68 2 : else serial=true;
69 418 : if(serial)log.printf(" doing calculation in serial\n");
70 418 : if( keywords.exists("LOWMEM") ) {
71 389 : plumed_assert( !keywords.exists("HIGHMEM") );
72 389 : parseFlag("LOWMEM",lowmem);
73 389 : if(lowmem) {
74 26 : log.printf(" lowering memory requirements\n");
75 26 : dertime_can_be_off=true;
76 : }
77 : }
78 418 : if( keywords.exists("HIGHMEM") ) {
79 21 : plumed_assert( !keywords.exists("LOWMEM") );
80 21 : bool highmem; parseFlag("HIGHMEM",highmem);
81 21 : lowmem=!highmem;
82 21 : if(!lowmem) log.printf(" increasing the memory requirements\n");
83 : }
84 418 : tolerance=nl_tolerance=epsilon;
85 418 : if( keywords.exists("TOL") ) parse("TOL",tolerance);
86 418 : if( tolerance>epsilon) {
87 6 : log.printf(" Ignoring contributions less than %f \n",tolerance);
88 : }
89 418 : parseFlag("TIMINGS",timers);
90 418 : stopwatch.start(); stopwatch.pause();
91 418 : }
92 :
93 836 : ActionWithVessel::~ActionWithVessel() {
94 418 : for(unsigned i=0; i<functions.size(); ++i) delete functions[i];
95 418 : stopwatch.start(); stopwatch.stop();
96 418 : 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 418 : delete &stopwatch;
101 418 : }
102 :
103 328 : void ActionWithVessel::addVessel( const std::string& name, const std::string& input, const int numlab ) {
104 328 : VesselOptions da(name,"",numlab,input,this);
105 328 : 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 484 : void ActionWithVessel::addVessel( Vessel* vv ) {
115 484 : ShortcutVessel* sv=dynamic_cast<ShortcutVessel*>(vv);
116 484 : if(!sv) { vv->checkRead(); functions.push_back(vv); }
117 494 : else { delete sv; return; }
118 :
119 474 : StoreDataVessel* mm=dynamic_cast<StoreDataVessel*>( vv );
120 474 : if( mydata && mm ) error("cannot have more than one StoreDataVessel in one action");
121 474 : else if( mm ) mydata=mm;
122 355 : else dertime_can_be_off=false;
123 : }
124 :
125 46 : BridgeVessel* ActionWithVessel::addBridgingVessel( ActionWithVessel* tome ) {
126 46 : 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 46 : 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 90 : 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 12645172 : void ActionWithVessel::addTaskToList( const unsigned& taskCode ) {
153 12645172 : fullTaskList.push_back( taskCode ); taskFlags.push_back(0);
154 12645172 : plumed_assert( fullTaskList.size()==taskFlags.size() );
155 12645172 : }
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 8620 : 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 8229 : if( vesselRegister().check(thiskey) ) {
166 2501 : 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 2408 : 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 22 : input.clear();
179 44 : }
180 : }
181 : }
182 8229 : }
183 :
184 : // Make sure all vessels have had been resized at start
185 391 : if( functions.size()>0 ) resizeFunctions();
186 391 : }
187 :
188 1270 : void ActionWithVessel::resizeFunctions() {
189 1270 : for(unsigned i=0; i<functions.size(); ++i) functions[i]->resize();
190 1270 : }
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 1046 : 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 3500 : void ActionWithVessel::lockContributors() {
205 3500 : nactive_tasks = 0;
206 18472065 : for(unsigned i=0; i<fullTaskList.size(); ++i) {
207 18468565 : if( taskFlags[i]>0 ) nactive_tasks++;
208 : }
209 :
210 3500 : unsigned n=0;
211 3500 : partialTaskList.resize( nactive_tasks );
212 3500 : indexOfTaskInFullList.resize( nactive_tasks );
213 18472065 : for(unsigned i=0; i<fullTaskList.size(); ++i) {
214 : // Deactivate sets inactive tasks to number not equal to zero
215 18468565 : if( taskFlags[i]>0 ) {
216 5509931 : partialTaskList[n] = fullTaskList[i];
217 5509931 : indexOfTaskInFullList[n]=i;
218 5509931 : n++;
219 : }
220 : }
221 : plumed_dbg_assert( n==nactive_tasks );
222 8090 : for(unsigned i=0; i<functions.size(); ++i) {
223 4590 : BridgeVessel* bb = dynamic_cast<BridgeVessel*>( functions[i] );
224 4590 : if( bb ) bb->copyTaskFlags();
225 : }
226 : // Resize mydata to accomodate all active tasks
227 3500 : if( mydata ) mydata->resize();
228 3500 : contributorsAreUnlocked=false;
229 3500 : }
230 :
231 3500 : void ActionWithVessel::deactivateAllTasks() {
232 3500 : contributorsAreUnlocked=true; nactive_tasks = 0;
233 3500 : taskFlags.assign(taskFlags.size(),0);
234 3500 : }
235 :
236 205651 : bool ActionWithVessel::taskIsCurrentlyActive( const unsigned& index ) const {
237 205651 : plumed_dbg_assert( index<taskFlags.size() ); return (taskFlags[index]>0);
238 : }
239 :
240 22875 : void ActionWithVessel::doJobsRequiredBeforeTaskList() {
241 : // Do any preparatory stuff for functions
242 22875 : for(unsigned j=0; j<functions.size(); ++j) functions[j]->prepare();
243 22875 : }
244 :
245 23388 : unsigned ActionWithVessel::getSizeOfBuffer( unsigned& bufsize ) {
246 23388 : for(unsigned i=0; i<functions.size(); ++i) functions[i]->setBufferStart( bufsize );
247 23388 : if( buffer.size()!=bufsize ) buffer.resize( bufsize );
248 23388 : if( mydata ) {
249 2315 : unsigned dsize=mydata->getSizeOfDerivativeList();
250 2315 : if( der_list.size()!=dsize ) der_list.resize( dsize );
251 : }
252 23388 : return bufsize;
253 : }
254 :
255 19850 : void ActionWithVessel::runAllTasks() {
256 19850 : plumed_massert( !contributorsAreUnlocked && functions.size()>0, "you must have a call to readVesselKeywords somewhere" );
257 19850 : unsigned stride=comm.Get_size();
258 19850 : unsigned rank=comm.Get_rank();
259 19850 : if(serial) { stride=1; rank=0; }
260 :
261 : // Make sure jobs are done
262 19850 : if(timers) stopwatch.start("1 Prepare Tasks");
263 19850 : doJobsRequiredBeforeTaskList();
264 19850 : if(timers) stopwatch.stop("1 Prepare Tasks");
265 :
266 : // Get number of threads for OpenMP
267 19850 : unsigned nt=OpenMP::getNumThreads();
268 19850 : if( nt*stride*2>nactive_tasks || !threadSafe()) nt=1;
269 :
270 : // Get size for buffer
271 19850 : unsigned bsize=0, bufsize=getSizeOfBuffer( bsize );
272 : // Clear buffer
273 19850 : buffer.assign( buffer.size(), 0.0 );
274 : // Switch off calculation of derivatives in main loop
275 19850 : if( dertime_can_be_off ) dertime=false;
276 :
277 19850 : if(timers) stopwatch.start("2 Loop over tasks");
278 48413 : #pragma omp parallel num_threads(nt)
279 : {
280 28563 : std::vector<double> omp_buffer;
281 47441 : if( nt>1 ) omp_buffer.resize( bufsize, 0.0 );
282 58371 : MultiValue myvals( getNumberOfQuantities(), getNumberOfDerivatives() );
283 58171 : MultiValue bvals( getNumberOfQuantities(), getNumberOfDerivatives() );
284 29340 : myvals.clearAll(); bvals.clearAll();
285 :
286 29275 : #pragma omp for nowait schedule(dynamic)
287 : for(unsigned i=rank; i<nactive_tasks; i+=stride) {
288 : // Calculate the stuff in the loop for this action
289 984713 : 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 984685 : if( myvals.get(0)<tolerance ) {
295 : // Clear the derivatives
296 577574 : myvals.clearAll();
297 577700 : 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 406708 : if( nt>1 ) {
304 363212 : calculateAllVessels( indexOfTaskInFullList[i], myvals, bvals, omp_buffer, der_list );
305 : } else {
306 43496 : calculateAllVessels( indexOfTaskInFullList[i], myvals, bvals, buffer, der_list );
307 : }
308 :
309 : // Clear the value
310 407010 : myvals.clearAll();
311 : }
312 58700 : #pragma omp critical
313 58690 : if(nt>1) for(unsigned i=0; i<bufsize; ++i) buffer[i]+=omp_buffer[i];
314 : }
315 19850 : if(timers) stopwatch.stop("2 Loop over tasks");
316 : // Turn back on derivative calculation
317 19850 : dertime=true;
318 :
319 19850 : if(timers) stopwatch.start("3 MPI gather");
320 : // MPI Gather everything
321 19850 : if( !serial && buffer.size()>0 ) comm.Sum( buffer );
322 : // MPI Gather index stores
323 19850 : 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 19850 : if(timers) stopwatch.stop("3 MPI gather");
329 :
330 19850 : if(timers) stopwatch.start("4 Finishing computations");
331 19850 : finishComputations( buffer );
332 19850 : if(timers) stopwatch.stop("4 Finishing computations");
333 19850 : }
334 :
335 0 : void ActionWithVessel::transformBridgedDerivatives( const unsigned& current, MultiValue& invals, MultiValue& outvals ) const {
336 0 : plumed_error();
337 : }
338 :
339 468998 : void ActionWithVessel::calculateAllVessels( const unsigned& taskCode, MultiValue& myvals, MultiValue& bvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) {
340 1070802 : 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 601539 : functions[j]->calculate( taskCode, functions[j]->transformDerivatives(taskCode, myvals, bvals), buffer, der_list );
344 601784 : if( !actionIsBridged ) bvals.clearAll();
345 : }
346 469502 : return;
347 : }
348 :
349 22875 : void ActionWithVessel::finishComputations( const std::vector<double>& buffer ) {
350 : // Set the final value of the function
351 22875 : for(unsigned j=0; j<functions.size(); ++j) functions[j]->finish( buffer );
352 22875 : }
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 6971 : forcesToApply.assign( forcesToApply.size(),0.0 );
361 6971 : bool wasforced=false;
362 20136 : for(unsigned i=0; i<getNumberOfVessels(); ++i) {
363 13165 : if( (functions[i]->applyForce( tmpforces )) ) {
364 278 : wasforced=true;
365 278 : 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 0 : 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 4821 : }
|
