Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2020 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 669 : void ActionWithVessel::registerKeywords(Keywords& keys) {
39 2676 : 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 in between 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 2676 : 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 2007 : keys.addFlag("SERIAL",false,"do the calculation in serial. Do not use MPI");
46 2007 : keys.addFlag("LOWMEM",false,"lower the memory requirements");
47 2007 : keys.addFlag("TIMINGS",false,"output information on the timings of the various parts of the calculation");
48 2007 : keys.reserveFlag("HIGHMEM",false,"use a more memory intensive version of this collective variable");
49 669 : keys.add( vesselRegister().getKeywords() );
50 669 : }
51 :
52 557 : 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 : dertime_can_be_off(false),
60 : dertime(true),
61 : contributorsAreUnlocked(false),
62 : weightHasDerivatives(false),
63 2228 : mydata(NULL)
64 : {
65 1114 : maxderivatives=309; parse("MAXDERIVATIVES",maxderivatives);
66 1639 : if( keywords.exists("SERIAL") ) parseFlag("SERIAL",serial);
67 32 : else serial=true;
68 557 : if(serial)log.printf(" doing calculation in serial\n");
69 1114 : if( keywords.exists("LOWMEM") ) {
70 1052 : plumed_assert( !keywords.exists("HIGHMEM") );
71 1052 : parseFlag("LOWMEM",lowmem);
72 526 : if(lowmem) {
73 29 : log.printf(" lowering memory requirements\n");
74 29 : dertime_can_be_off=true;
75 : }
76 : }
77 1114 : if( keywords.exists("HIGHMEM") ) {
78 46 : plumed_assert( !keywords.exists("LOWMEM") );
79 46 : bool highmem; parseFlag("HIGHMEM",highmem);
80 23 : lowmem=!highmem;
81 23 : if(!lowmem) log.printf(" increasing the memory requirements\n");
82 : }
83 557 : tolerance=nl_tolerance=epsilon;
84 1560 : if( keywords.exists("TOL") ) parse("TOL",tolerance);
85 557 : if( tolerance>epsilon) {
86 6 : log.printf(" Ignoring contributions less than %f \n",tolerance);
87 : }
88 1114 : parseFlag("TIMINGS",timers);
89 557 : stopwatch.start(); stopwatch.pause();
90 557 : }
91 :
92 1114 : ActionWithVessel::~ActionWithVessel() {
93 557 : stopwatch.start(); stopwatch.stop();
94 557 : if(timers) {
95 0 : log.printf("timings for action %s with label %s \n", getName().c_str(), getLabel().c_str() );
96 0 : log<<stopwatch;
97 : }
98 557 : }
99 :
100 342 : void ActionWithVessel::addVessel( const std::string& name, const std::string& input, const int numlab ) {
101 1026 : VesselOptions da(name,"",numlab,input,this);
102 1026 : auto vv=vesselRegister().create(name,da);
103 342 : FunctionVessel* fv=dynamic_cast<FunctionVessel*>(vv.get());
104 342 : if( fv ) {
105 312 : std::string mylabel=Vessel::transformName( name );
106 312 : plumed_massert( keywords.outputComponentExists(mylabel,false), "a description of the value calculated by vessel " + name + " has not been added to the manual");
107 : }
108 684 : addVessel(std::move(vv));
109 342 : }
110 :
111 535 : void ActionWithVessel::addVessel( std::unique_ptr<Vessel> vv_ptr ) {
112 :
113 : // In the original code, the dynamically casted pointer was deleted here.
114 : // Now that vv_ptr is a unique_ptr, the object will be deleted automatically when
115 : // exiting this routine.
116 535 : if(dynamic_cast<ShortcutVessel*>(vv_ptr.get())) return;
117 :
118 524 : vv_ptr->checkRead();
119 :
120 524 : StoreDataVessel* mm=dynamic_cast<StoreDataVessel*>( vv_ptr.get() );
121 524 : if( mydata && mm ) error("cannot have more than one StoreDataVessel in one action");
122 524 : else if( mm ) mydata=mm;
123 393 : else dertime_can_be_off=false;
124 :
125 : // Ownership is transfered to functions
126 524 : functions.emplace_back(std::move(vv_ptr));
127 : }
128 :
129 46 : BridgeVessel* ActionWithVessel::addBridgingVessel( ActionWithVessel* tome ) {
130 230 : VesselOptions da("","",0,"",this);
131 46 : std::unique_ptr<BridgeVessel> bv(new BridgeVessel(da));
132 46 : bv->setOutputAction( tome );
133 46 : tome->actionIsBridged=true; dertime_can_be_off=false;
134 : // store this pointer in order to return it later.
135 : // notice that I cannot access this with functions.tail().get()
136 : // since functions contains pointers to a different class (Vessel)
137 : auto toBeReturned=bv.get();
138 46 : functions.emplace_back( std::move(bv) );
139 46 : resizeFunctions();
140 46 : return toBeReturned;
141 : }
142 :
143 192 : StoreDataVessel* ActionWithVessel::buildDataStashes( ActionWithVessel* actionThatUses ) {
144 192 : if(mydata) {
145 87 : if( actionThatUses ) mydata->addActionThatUses( actionThatUses );
146 87 : return mydata;
147 : }
148 :
149 525 : VesselOptions da("","",0,"",this);
150 105 : std::unique_ptr<StoreDataVessel> mm( new StoreDataVessel(da) );
151 105 : if( actionThatUses ) mm->addActionThatUses( actionThatUses );
152 210 : addVessel(std::move(mm));
153 :
154 : // Make sure resizing of vessels is done
155 105 : resizeFunctions();
156 :
157 105 : return mydata;
158 : }
159 :
160 12657137 : void ActionWithVessel::addTaskToList( const unsigned& taskCode ) {
161 25314274 : fullTaskList.push_back( taskCode ); taskFlags.push_back(0);
162 12657137 : plumed_assert( fullTaskList.size()==taskFlags.size() );
163 12657137 : }
164 :
165 407 : void ActionWithVessel::readVesselKeywords() {
166 : // Set maxderivatives if it is too big
167 407 : if( maxderivatives>getNumberOfDerivatives() ) maxderivatives=getNumberOfDerivatives();
168 :
169 : // Loop over all keywords find the vessels and create appropriate functions
170 18391 : for(unsigned i=0; i<keywords.size(); ++i) {
171 17984 : std::string thiskey,input; thiskey=keywords.getKeyword(i);
172 : // Check if this is a key for a vessel
173 26976 : if( vesselRegister().check(thiskey) ) {
174 5306 : plumed_assert( keywords.style(thiskey,"vessel") );
175 2653 : bool dothis=false; parseFlag(thiskey,dothis);
176 2653 : if(dothis) addVessel( thiskey, input );
177 :
178 2653 : parse(thiskey,input);
179 2653 : if(input.size()!=0) {
180 118 : addVessel( thiskey, input );
181 : } else {
182 25 : for(unsigned i=1;; ++i) {
183 2560 : if( !parseNumbered(thiskey,i,input) ) break;
184 25 : std::string ss; Tools::convert(i,ss);
185 25 : addVessel( thiskey, input, i );
186 : input.clear();
187 25 : }
188 : }
189 : }
190 : }
191 :
192 : // Make sure all vessels have had been resized at start
193 407 : if( functions.size()>0 ) resizeFunctions();
194 407 : }
195 :
196 1359 : void ActionWithVessel::resizeFunctions() {
197 9168 : for(unsigned i=0; i<functions.size(); ++i) functions[i]->resize();
198 1359 : }
199 :
200 809 : void ActionWithVessel::needsDerivatives() {
201 : // Turn on the derivatives and resize
202 809 : noderiv=false; resizeFunctions();
203 : // Setting contributors unlocked here ensures that link cells are ignored
204 809 : contributorsAreUnlocked=true; contributorsAreUnlocked=false;
205 : // And turn on the derivatives in all actions on which we are dependent
206 2512 : for(unsigned i=0; i<getDependencies().size(); ++i) {
207 298 : ActionWithVessel* vv=dynamic_cast<ActionWithVessel*>( getDependencies()[i] );
208 298 : if(vv) vv->needsDerivatives();
209 : }
210 809 : }
211 :
212 3302 : void ActionWithVessel::lockContributors() {
213 3302 : nactive_tasks = 0;
214 55462303 : for(unsigned i=0; i<fullTaskList.size(); ++i) {
215 18485233 : if( taskFlags[i]>0 ) nactive_tasks++;
216 : }
217 :
218 : unsigned n=0;
219 3302 : partialTaskList.resize( nactive_tasks );
220 3302 : indexOfTaskInFullList.resize( nactive_tasks );
221 55462303 : for(unsigned i=0; i<fullTaskList.size(); ++i) {
222 : // Deactivate sets inactive tasks to number not equal to zero
223 18485233 : if( taskFlags[i]>0 ) {
224 10956464 : partialTaskList[n] = fullTaskList[i];
225 5478232 : indexOfTaskInFullList[n]=i;
226 5478232 : n++;
227 : }
228 : }
229 : plumed_dbg_assert( n==nactive_tasks );
230 21181 : for(unsigned i=0; i<functions.size(); ++i) {
231 4859 : BridgeVessel* bb = dynamic_cast<BridgeVessel*>( functions[i].get() );
232 4859 : if( bb ) bb->copyTaskFlags();
233 : }
234 : // Resize mydata to accomodate all active tasks
235 3302 : if( mydata ) mydata->resize();
236 3302 : contributorsAreUnlocked=false;
237 3302 : }
238 :
239 3302 : void ActionWithVessel::deactivateAllTasks() {
240 3302 : contributorsAreUnlocked=true; nactive_tasks = 0;
241 6604 : taskFlags.assign(taskFlags.size(),0);
242 3302 : }
243 :
244 213721 : bool ActionWithVessel::taskIsCurrentlyActive( const unsigned& index ) const {
245 427442 : plumed_dbg_assert( index<taskFlags.size() ); return (taskFlags[index]>0);
246 : }
247 :
248 22680 : void ActionWithVessel::doJobsRequiredBeforeTaskList() {
249 : // Do any preparatory stuff for functions
250 150213 : for(unsigned j=0; j<functions.size(); ++j) functions[j]->prepare();
251 22680 : }
252 :
253 23193 : unsigned ActionWithVessel::getSizeOfBuffer( unsigned& bufsize ) {
254 153267 : for(unsigned i=0; i<functions.size(); ++i) functions[i]->setBufferStart( bufsize );
255 23193 : if( buffer.size()!=bufsize ) buffer.resize( bufsize );
256 23193 : if( mydata ) {
257 : unsigned dsize=mydata->getSizeOfDerivativeList();
258 2553 : if( der_list.size()!=dsize ) der_list.resize( dsize );
259 : }
260 23193 : return bufsize;
261 : }
262 :
263 19655 : void ActionWithVessel::runAllTasks() {
264 39310 : plumed_massert( !contributorsAreUnlocked && functions.size()>0, "you must have a call to readVesselKeywords somewhere" );
265 19655 : unsigned stride=comm.Get_size();
266 19655 : unsigned rank=comm.Get_rank();
267 19655 : if(serial) { stride=1; rank=0; }
268 :
269 : // Make sure jobs are done
270 19655 : if(timers) stopwatch.start("1 Prepare Tasks");
271 19655 : doJobsRequiredBeforeTaskList();
272 19655 : if(timers) stopwatch.stop("1 Prepare Tasks");
273 :
274 : // Get number of threads for OpenMP
275 19655 : unsigned nt=OpenMP::getNumThreads();
276 19655 : if( nt*stride*2>nactive_tasks || !threadSafe()) nt=1;
277 :
278 : // Get size for buffer
279 19655 : unsigned bsize=0, bufsize=getSizeOfBuffer( bsize );
280 : // Clear buffer
281 39310 : buffer.assign( buffer.size(), 0.0 );
282 : // Switch off calculation of derivatives in main loop
283 19655 : if( dertime_can_be_off ) dertime=false;
284 :
285 19655 : if(timers) stopwatch.start("2 Loop over tasks");
286 48574 : #pragma omp parallel num_threads(nt)
287 : {
288 : std::vector<double> omp_buffer;
289 47447 : if( nt>1 ) omp_buffer.resize( bufsize, 0.0 );
290 57838 : MultiValue myvals( getNumberOfQuantities(), getNumberOfDerivatives() );
291 57838 : MultiValue bvals( getNumberOfQuantities(), getNumberOfDerivatives() );
292 28919 : myvals.clearAll(); bvals.clearAll();
293 :
294 : #pragma omp for nowait schedule(dynamic)
295 28919 : for(unsigned i=rank; i<nactive_tasks; i+=stride) {
296 : // Calculate the stuff in the loop for this action
297 1899540 : performTask( indexOfTaskInFullList[i], partialTaskList[i], myvals );
298 :
299 : // Check for conditions that allow us to just to skip the calculation
300 : // the condition is that the weight of the contribution is low
301 : // N.B. Here weights are assumed to be between zero and one
302 1527608 : if( myvals.get(0)<tolerance ) {
303 : // Clear the derivatives
304 577838 : myvals.clearAll();
305 577838 : continue;
306 : }
307 :
308 : // Now calculate all the functions
309 : // If the contribution of this quantity is very small at neighbour list time ignore it
310 : // untill next neighbour list time
311 371932 : if( nt>1 ) {
312 656464 : calculateAllVessels( indexOfTaskInFullList[i], myvals, bvals, omp_buffer, der_list );
313 : } else {
314 87400 : calculateAllVessels( indexOfTaskInFullList[i], myvals, bvals, buffer, der_list );
315 : }
316 :
317 : // Clear the value
318 371932 : myvals.clearAll();
319 : }
320 57838 : #pragma omp critical
321 38747735 : if(nt>1) for(unsigned i=0; i<bufsize; ++i) buffer[i]+=omp_buffer[i];
322 : }
323 19655 : if(timers) stopwatch.stop("2 Loop over tasks");
324 : // Turn back on derivative calculation
325 19655 : dertime=true;
326 :
327 19655 : if(timers) stopwatch.start("3 MPI gather");
328 : // MPI Gather everything
329 39310 : if( !serial && buffer.size()>0 ) comm.Sum( buffer );
330 : // MPI Gather index stores
331 19655 : if( mydata && !lowmem && !noderiv ) {
332 689 : comm.Sum( der_list ); mydata->setActiveValsAndDerivatives( der_list );
333 : }
334 : // Update the elements that are makign contributions to the sum here
335 : // this causes problems if we do it in prepare
336 19655 : if(timers) stopwatch.stop("3 MPI gather");
337 :
338 19655 : if(timers) stopwatch.start("4 Finishing computations");
339 19655 : finishComputations( buffer );
340 19655 : if(timers) stopwatch.stop("4 Finishing computations");
341 19655 : }
342 :
343 0 : void ActionWithVessel::transformBridgedDerivatives( const unsigned& current, MultiValue& invals, MultiValue& outvals ) const {
344 0 : plumed_error();
345 : }
346 :
347 434440 : void ActionWithVessel::calculateAllVessels( const unsigned& taskCode, MultiValue& myvals, MultiValue& bvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) {
348 2602754 : for(unsigned j=0; j<functions.size(); ++j) {
349 : // Calculate returns a bool that tells us if this particular
350 : // quantity is contributing more than the tolerance
351 577958 : functions[j]->calculate( taskCode, functions[j]->transformDerivatives(taskCode, myvals, bvals), buffer, der_list );
352 577958 : if( !actionIsBridged ) bvals.clearAll();
353 : }
354 434440 : return;
355 : }
356 :
357 22680 : void ActionWithVessel::finishComputations( const std::vector<double>& buffer ) {
358 : // Set the final value of the function
359 150213 : for(unsigned j=0; j<functions.size(); ++j) functions[j]->finish( buffer );
360 22680 : }
361 :
362 7210 : bool ActionWithVessel::getForcesFromVessels( std::vector<double>& forcesToApply ) {
363 : #ifndef NDEBUG
364 : if( forcesToApply.size()>0 ) plumed_dbg_assert( forcesToApply.size()==getNumberOfDerivatives() );
365 : #endif
366 7210 : if(tmpforces.size()!=forcesToApply.size() ) tmpforces.resize( forcesToApply.size() );
367 :
368 14420 : forcesToApply.assign( forcesToApply.size(),0.0 );
369 : bool wasforced=false;
370 34948 : for(unsigned i=0; i<getNumberOfVessels(); ++i) {
371 27738 : if( (functions[i]->applyForce( tmpforces )) ) {
372 : wasforced=true;
373 2096364 : for(unsigned j=0; j<forcesToApply.size(); ++j) forcesToApply[j]+=tmpforces[j];
374 : }
375 : }
376 7210 : return wasforced;
377 : }
378 :
379 0 : void ActionWithVessel::retrieveDomain( std::string& min, std::string& max ) {
380 0 : plumed_merror("If your function is periodic you need to add a retrieveDomain function so that ActionWithVessel can retrieve the domain");
381 : }
382 :
383 0 : Vessel* ActionWithVessel::getVesselWithName( const std::string& mynam ) {
384 : int target=-1;
385 0 : for(unsigned i=0; i<functions.size(); ++i) {
386 0 : if( functions[i]->getName().find(mynam)!=std::string::npos ) {
387 0 : if( target<0 ) target=i;
388 0 : else error("found more than one " + mynam + " object in action");
389 : }
390 : }
391 0 : return functions[target].get();
392 : }
393 :
394 : }
395 5517 : }
|
