Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2018 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 414 : 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 414 : "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 414 : "we have to start using lowmem");
45 414 : keys.addFlag("SERIAL",false,"do the calculation in serial. Do not parallelize");
46 414 : keys.addFlag("LOWMEM",false,"lower the memory requirements");
47 414 : keys.addFlag("TIMINGS",false,"output information on the timings of the various parts of the calculation");
48 414 : keys.reserveFlag("HIGHMEM",false,"use a more memory intensive version of this collective variable");
49 414 : keys.add( vesselRegister().getKeywords() );
50 414 : }
51 :
52 341 : 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 341 : stopwatch(*new Stopwatch),
60 : dertime_can_be_off(false),
61 : dertime(true),
62 : contributorsAreUnlocked(false),
63 : weightHasDerivatives(false),
64 682 : mydata(NULL)
65 : {
66 341 : maxderivatives=309; parse("MAXDERIVATIVES",maxderivatives);
67 341 : if( keywords.exists("SERIAL") ) parseFlag("SERIAL",serial);
68 2 : else serial=true;
69 341 : if(serial)log.printf(" doing calculation in serial\n");
70 341 : if( keywords.exists("LOWMEM") ) {
71 320 : plumed_assert( !keywords.exists("HIGHMEM") );
72 320 : parseFlag("LOWMEM",lowmem);
73 320 : if(lowmem) {
74 16 : log.printf(" lowering memory requirements\n");
75 16 : dertime_can_be_off=true;
76 : }
77 : }
78 341 : if( keywords.exists("HIGHMEM") ) {
79 15 : plumed_assert( !keywords.exists("LOWMEM") );
80 15 : bool highmem; parseFlag("HIGHMEM",highmem);
81 15 : lowmem=!highmem;
82 15 : if(!lowmem) log.printf(" increasing the memory requirements\n");
83 : }
84 341 : tolerance=nl_tolerance=epsilon;
85 341 : if( keywords.exists("TOL") ) parse("TOL",tolerance);
86 341 : if( tolerance>epsilon) {
87 6 : log.printf(" Ignoring contributions less than %f \n",tolerance);
88 : }
89 341 : parseFlag("TIMINGS",timers);
90 341 : stopwatch.start(); stopwatch.pause();
91 341 : }
92 :
93 682 : ActionWithVessel::~ActionWithVessel() {
94 341 : for(unsigned i=0; i<functions.size(); ++i) delete functions[i];
95 341 : stopwatch.start(); stopwatch.stop();
96 341 : 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 341 : delete &stopwatch;
101 341 : }
102 :
103 301 : void ActionWithVessel::addVessel( const std::string& name, const std::string& input, const int numlab ) {
104 301 : VesselOptions da(name,"",numlab,input,this);
105 301 : Vessel* vv=vesselRegister().create(name,da);
106 301 : FunctionVessel* fv=dynamic_cast<FunctionVessel*>(vv);
107 301 : if( fv ) {
108 279 : std::string mylabel=Vessel::transformName( name );
109 279 : plumed_massert( keywords.outputComponentExists(mylabel,false), "a description of the value calculated by vessel " + name + " has not been added to the manual");
110 : }
111 301 : addVessel(vv);
112 301 : }
113 :
114 423 : void ActionWithVessel::addVessel( Vessel* vv ) {
115 423 : ShortcutVessel* sv=dynamic_cast<ShortcutVessel*>(vv);
116 423 : if(!sv) { vv->checkRead(); functions.push_back(vv); }
117 433 : else { delete sv; return; }
118 :
119 413 : StoreDataVessel* mm=dynamic_cast<StoreDataVessel*>( vv );
120 413 : if( mydata && mm ) error("cannot have more than one StoreDataVessel in one action");
121 413 : else if( mm ) mydata=mm;
122 324 : else dertime_can_be_off=false;
123 : }
124 :
125 45 : BridgeVessel* ActionWithVessel::addBridgingVessel( ActionWithVessel* tome ) {
126 45 : VesselOptions da("","",0,"",this);
127 45 : BridgeVessel* bv=new BridgeVessel(da);
128 45 : bv->setOutputAction( tome );
129 45 : tome->actionIsBridged=true; dertime_can_be_off=false;
130 45 : functions.push_back( dynamic_cast<Vessel*>(bv) );
131 45 : resizeFunctions();
132 45 : return bv;
133 : }
134 :
135 122 : StoreDataVessel* ActionWithVessel::buildDataStashes( ActionWithVessel* actionThatUses ) {
136 122 : if(mydata) {
137 53 : if( actionThatUses ) mydata->addActionThatUses( actionThatUses );
138 53 : return mydata;
139 : }
140 :
141 69 : VesselOptions da("","",0,"",this);
142 69 : StoreDataVessel* mm=new StoreDataVessel(da);
143 69 : if( actionThatUses ) mm->addActionThatUses( actionThatUses );
144 69 : addVessel(mm);
145 :
146 : // Make sure resizing of vessels is done
147 69 : resizeFunctions();
148 :
149 69 : return mydata;
150 : }
151 :
152 12401211 : void ActionWithVessel::addTaskToList( const unsigned& taskCode ) {
153 12401211 : fullTaskList.push_back( taskCode ); taskFlags.push_back(0);
154 12401211 : plumed_assert( fullTaskList.size()==taskFlags.size() );
155 12401211 : }
156 :
157 321 : void ActionWithVessel::readVesselKeywords() {
158 : // Set maxderivatives if it is too big
159 321 : if( maxderivatives>getNumberOfDerivatives() ) maxderivatives=getNumberOfDerivatives();
160 :
161 : // Loop over all keywords find the vessels and create appropriate functions
162 7119 : for(unsigned i=0; i<keywords.size(); ++i) {
163 13596 : std::string thiskey,input; thiskey=keywords.getKeyword(i);
164 : // Check if this is a key for a vessel
165 6798 : if( vesselRegister().check(thiskey) ) {
166 2119 : plumed_assert( keywords.style(thiskey,"vessel") );
167 2119 : bool dothis=false; parseFlag(thiskey,dothis);
168 2119 : if(dothis) addVessel( thiskey, input );
169 :
170 2119 : parse(thiskey,input);
171 2119 : if(input.size()!=0) {
172 112 : addVessel( thiskey, input );
173 : } else {
174 2029 : for(unsigned i=1;; ++i) {
175 2029 : 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 6798 : }
183 :
184 : // Make sure all vessels have had been resized at start
185 321 : if( functions.size()>0 ) resizeFunctions();
186 321 : }
187 :
188 1041 : void ActionWithVessel::resizeFunctions() {
189 1041 : for(unsigned i=0; i<functions.size(); ++i) functions[i]->resize();
190 1041 : }
191 :
192 613 : void ActionWithVessel::needsDerivatives() {
193 : // Turn on the derivatives and resize
194 613 : noderiv=false; resizeFunctions();
195 : // Setting contributors unlocked here ensures that link cells are ignored
196 613 : contributorsAreUnlocked=true; contributorsAreUnlocked=false;
197 : // And turn on the derivatives in all actions on which we are dependent
198 826 : for(unsigned i=0; i<getDependencies().size(); ++i) {
199 213 : ActionWithVessel* vv=dynamic_cast<ActionWithVessel*>( getDependencies()[i] );
200 213 : if(vv) vv->needsDerivatives();
201 : }
202 613 : }
203 :
204 3208 : void ActionWithVessel::lockContributors() {
205 3208 : nactive_tasks = 0;
206 17759341 : for(unsigned i=0; i<fullTaskList.size(); ++i) {
207 17756133 : if( taskFlags[i]>0 ) nactive_tasks++;
208 : }
209 :
210 3208 : unsigned n=0;
211 3208 : partialTaskList.resize( nactive_tasks );
212 3208 : indexOfTaskInFullList.resize( nactive_tasks );
213 17759341 : for(unsigned i=0; i<fullTaskList.size(); ++i) {
214 : // Deactivate sets inactive tasks to number not equal to zero
215 17756133 : if( taskFlags[i]>0 ) {
216 5150655 : partialTaskList[n] = fullTaskList[i];
217 5150655 : indexOfTaskInFullList[n]=i;
218 5150655 : n++;
219 : }
220 : }
221 : plumed_dbg_assert( n==nactive_tasks );
222 7488 : for(unsigned i=0; i<functions.size(); ++i) {
223 4280 : BridgeVessel* bb = dynamic_cast<BridgeVessel*>( functions[i] );
224 4280 : if( bb ) bb->copyTaskFlags();
225 : }
226 : // Resize mydata to accomodate all active tasks
227 3208 : if( mydata ) mydata->resize();
228 3208 : contributorsAreUnlocked=false;
229 3208 : }
230 :
231 3208 : void ActionWithVessel::deactivateAllTasks() {
232 3208 : contributorsAreUnlocked=true; nactive_tasks = 0;
233 3208 : taskFlags.assign(taskFlags.size(),0);
234 3208 : }
235 :
236 111680 : bool ActionWithVessel::taskIsCurrentlyActive( const unsigned& index ) const {
237 111680 : plumed_dbg_assert( index<taskFlags.size() ); return (taskFlags[index]>0);
238 : }
239 :
240 19548 : void ActionWithVessel::doJobsRequiredBeforeTaskList() {
241 : // Do any preparatory stuff for functions
242 19548 : for(unsigned j=0; j<functions.size(); ++j) functions[j]->prepare();
243 19548 : }
244 :
245 20056 : unsigned ActionWithVessel::getSizeOfBuffer( unsigned& bufsize ) {
246 20056 : for(unsigned i=0; i<functions.size(); ++i) functions[i]->setBufferStart( bufsize );
247 20056 : if( buffer.size()!=bufsize ) buffer.resize( bufsize );
248 20056 : if( mydata ) {
249 1012 : unsigned dsize=mydata->getSizeOfDerivativeList();
250 1012 : if( der_list.size()!=dsize ) der_list.resize( dsize );
251 : }
252 20056 : return bufsize;
253 : }
254 :
255 16544 : void ActionWithVessel::runAllTasks() {
256 16544 : plumed_massert( !contributorsAreUnlocked && functions.size()>0, "you must have a call to readVesselKeywords somewhere" );
257 16544 : unsigned stride=comm.Get_size();
258 16544 : unsigned rank=comm.Get_rank();
259 16544 : if(serial) { stride=1; rank=0; }
260 :
261 : // Make sure jobs are done
262 16544 : if(timers) stopwatch.start("1 Prepare Tasks");
263 16544 : doJobsRequiredBeforeTaskList();
264 16544 : if(timers) stopwatch.stop("1 Prepare Tasks");
265 :
266 : // Get number of threads for OpenMP
267 16544 : unsigned nt=OpenMP::getNumThreads();
268 16544 : if( nt*stride*2>nactive_tasks || !threadSafe()) nt=1;
269 :
270 : // Get size for buffer
271 16544 : unsigned bsize=0, bufsize=getSizeOfBuffer( bsize );
272 : // Clear buffer
273 16544 : buffer.assign( buffer.size(), 0.0 );
274 : // Switch off calculation of derivatives in main loop
275 16544 : if( dertime_can_be_off ) dertime=false;
276 :
277 16544 : if(timers) stopwatch.start("2 Loop over tasks");
278 39701 : #pragma omp parallel num_threads(nt)
279 : {
280 23157 : std::vector<double> omp_buffer;
281 38130 : if( nt>1 ) omp_buffer.resize( bufsize, 0.0 );
282 49956 : MultiValue myvals( getNumberOfQuantities(), getNumberOfDerivatives() );
283 47806 : MultiValue bvals( getNumberOfQuantities(), getNumberOfDerivatives() );
284 24500 : myvals.clearAll(); bvals.clearAll();
285 :
286 24445 : #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 875626 : 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 874468 : if( myvals.get(0)<tolerance ) {
295 : // Clear the derivatives
296 555623 : myvals.clearAll();
297 555748 : 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 319258 : if( nt>1 ) {
304 278580 : calculateAllVessels( indexOfTaskInFullList[i], myvals, bvals, omp_buffer, der_list );
305 : } else {
306 40678 : calculateAllVessels( indexOfTaskInFullList[i], myvals, bvals, buffer, der_list );
307 : }
308 :
309 : // Clear the value
310 319301 : myvals.clearAll();
311 : }
312 47768 : #pragma omp critical
313 48392 : if(nt>1) for(unsigned i=0; i<bufsize; ++i) buffer[i]+=omp_buffer[i];
314 : }
315 16544 : if(timers) stopwatch.stop("2 Loop over tasks");
316 : // Turn back on derivative calculation
317 16544 : dertime=true;
318 :
319 16544 : if(timers) stopwatch.start("3 MPI gather");
320 : // MPI Gather everything
321 16544 : if( !serial && buffer.size()>0 ) comm.Sum( buffer );
322 : // MPI Gather index stores
323 16544 : if( mydata && !lowmem && !noderiv ) {
324 631 : 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 16544 : if(timers) stopwatch.stop("3 MPI gather");
329 :
330 16544 : if(timers) stopwatch.start("4 Finishing computations");
331 16544 : finishComputations( buffer );
332 16544 : if(timers) stopwatch.stop("4 Finishing computations");
333 16544 : }
334 :
335 0 : void ActionWithVessel::transformBridgedDerivatives( const unsigned& current, MultiValue& invals, MultiValue& outvals ) const {
336 0 : plumed_error();
337 : }
338 :
339 380972 : void ActionWithVessel::calculateAllVessels( const unsigned& taskCode, MultiValue& myvals, MultiValue& bvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) {
340 892303 : 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 511056 : functions[j]->calculate( taskCode, functions[j]->transformDerivatives(taskCode, myvals, bvals), buffer, der_list );
344 511252 : if( !actionIsBridged ) bvals.clearAll();
345 : }
346 381420 : return;
347 : }
348 :
349 19548 : void ActionWithVessel::finishComputations( const std::vector<double>& buffer ) {
350 : // Set the final value of the function
351 19548 : for(unsigned j=0; j<functions.size(); ++j) functions[j]->finish( buffer );
352 19548 : }
353 :
354 5550 : bool ActionWithVessel::getForcesFromVessels( std::vector<double>& forcesToApply ) {
355 : #ifndef NDEBUG
356 : if( forcesToApply.size()>0 ) plumed_dbg_assert( forcesToApply.size()==getNumberOfDerivatives() );
357 : #endif
358 5550 : if(tmpforces.size()!=forcesToApply.size() ) tmpforces.resize( forcesToApply.size() );
359 :
360 5550 : forcesToApply.assign( forcesToApply.size(),0.0 );
361 5550 : bool wasforced=false;
362 17293 : for(unsigned i=0; i<getNumberOfVessels(); ++i) {
363 11743 : if( (functions[i]->applyForce( tmpforces )) ) {
364 145 : wasforced=true;
365 145 : for(unsigned j=0; j<forcesToApply.size(); ++j) forcesToApply[j]+=tmpforces[j];
366 : }
367 : }
368 5550 : 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 2523 : }
|
