Loading [MathJax]/extensions/TeX/AMSsymbols.js
PSCF v1.2
rpg/fts/montecarlo/McSimulator.tpp
1#ifndef RPG_MC_SIMULATOR_TPP
2#define RPG_MC_SIMULATOR_TPP
3
4/*
5* PSCF - Polymer Self-Consistent Field Theory
6*
7* Copyright 2016 - 2022, The Regents of the University of Minnesota
8* Distributed under the terms of the GNU General Public License.
9*/
10
11#include "McSimulator.h"
12
13#include <rpg/System.h>
14#include <rpg/fts/montecarlo/McMoveFactory.h>
15#include <rpg/fts/analyzer/AnalyzerFactory.h>
16#include <rpg/fts/trajectory/TrajectoryReader.h>
17#include <rpg/fts/trajectory/TrajectoryReaderFactory.h>
18#include <rpg/fts/compressor/Compressor.h>
19#include <rpg/fts/perturbation/Perturbation.h>
20#include <rpg/fts/perturbation/PerturbationFactory.h>
21#include <rpg/fts/ramp/RampFactory.h>
22#include <rpg/fts/ramp/Ramp.h>
23
24#include <util/random/Random.h>
25#include <pscf/cuda/CudaRandom.h>
26#include <util/misc/Timer.h>
27#include <util/global.h>
28
29#include <gsl/gsl_eigen.h>
30
31namespace Pscf {
32namespace Rpg {
33
34 using namespace Util;
35 using namespace Pscf::Prdc::Cuda;
36
37 /*
38 * Constructor.
39 */
40 template <int D>
42 : Simulator<D>(system),
43 mcMoveManager_(*this, system),
44 analyzerManager_(*this, system),
45 trajectoryReaderFactoryPtr_(0)
46 {
47 setClassName("McSimulator");
48 trajectoryReaderFactoryPtr_
50 }
51
52 /*
53 * Destructor.
54 */
55 template <int D>
57 {
58 if (trajectoryReaderFactoryPtr_) {
59 delete trajectoryReaderFactoryPtr_;
60 }
61 }
62
63 /*
64 * Read instructions for creating objects from file.
65 */
66 template <int D>
67 void McSimulator<D>::readParameters(std::istream &in)
68 {
69 // Optionally read random seed, initialize random number generators
70 readRandomSeed(in);
71
72 // Attempt to read optional McMoveManager block
73 readParamCompositeOptional(in, mcMoveManager_);
74
75 // Attempt to read optional Compressor, Perturbation and Ramp blocks
76 bool isEnd = false;
77 readCompressor(in, isEnd);
78 readPerturbation(in, isEnd);
79 if (hasMcMoves()) {
80 UTIL_CHECK(hasCompressor());
81 readRamp(in, isEnd);
82 }
83 // Tthe Compressor is required if McMoves are declared
84 // Ramp is allowed only if McMoves are declared
85
86 // Optionally read AnalyzerManager block
87 Analyzer<D>::baseInterval = 0; // default value
88 readParamCompositeOptional(in, analyzerManager_);
89
90 // Figure out what needs to be saved in stored state
91 state_.needsCc = false;
92 state_.needsDc = false;
93 state_.needsHamiltonian = true;
94 if (mcMoveManager_.needsCc()){
95 state_.needsCc = true;
96 }
97 if (mcMoveManager_.needsDc()){
98 state_.needsDc = true;
99 }
100
101 // Initialize Simulator<D> base class
102 allocate();
103 }
104
105 /*
106 * Initialize just prior to a run.
107 */
108 template <int D>
109 void McSimulator<D>::setup(int nStep)
110 {
111 UTIL_CHECK(system().w().hasData());
112
113 // Eigenanalysis of the projected chi matrix.
114 analyzeChi();
115
116 if (hasPerturbation()) {
117 perturbation().setup();
118 }
119
120 if (hasRamp()) {
121 ramp().setup(nStep);
122 }
123
124 // Solve the MDE for the initial state
125 system().compute();
126
127 // Compress the initial field before entering simulate loop.
128 if (hasCompressor()) {
129 compressor().compress();
130 compressor().clearTimers();
131 }
132
133 // Compute field components and Hamiltonian
134 computeWc();
135 if (state_.needsCc || state_.needsDc) {
136 computeCc();
137 }
138 if (state_.needsDc) {
139 computeDc();
140 }
141 computeHamiltonian();
142
143 // Setup MC moves (if any)
144 if (hasMcMoves()) {
145 mcMoveManager_.setup();
146 }
147
148 // Setup analyzers (if any)
149 if (analyzerManager_.size() > 0){
150 analyzerManager_.setup();
151 }
152
153 }
154
155 /*
156 * Perform a field theoretic MC simulation of nStep steps.
157 */
158 template <int D>
160 {
161 UTIL_CHECK(hasMcMoves());
162 UTIL_CHECK(hasCompressor());
163
164 // Initial setup
165 setup(nStep);
166
167 Log::file() << std::endl;
168
169 // Main Monte Carlo loop
170 Timer timer;
171 Timer analyzerTimer;
172 timer.start();
173 iStep_ = 0;
174 if (hasRamp()) {
175 ramp().setParameters(iStep_);
176 }
177
178 // Analysis initial step (if any)
179 analyzerTimer.start();
180 analyzerManager_.sample(iStep_);
181 analyzerTimer.stop();
182
183 for (iTotalStep_ = 0; iTotalStep_ < nStep; ++iTotalStep_) {
184 // Choose and attempt an McMove
185 bool converged;
186 converged = mcMoveManager_.chooseMove().move();
187
188 if (converged){
189 iStep_++;
190
191 if (hasRamp()) {
192 ramp().setParameters(iStep_);
193 }
194
195 // Analysis (if any)
196 analyzerTimer.start();
197 if (Analyzer<D>::baseInterval != 0) {
198 if (iStep_ % Analyzer<D>::baseInterval == 0) {
199 if (analyzerManager_.size() > 0) {
200 analyzerManager_.sample(iStep_);
201 }
202 }
203 }
204 analyzerTimer.stop();
205
206 } else{
207 Log::file() << "Step: "<< iTotalStep_<< " fail to converge" << "\n";
208 }
209 }
210
211 timer.stop();
212 double time = timer.time();
213 double analyzerTime = analyzerTimer.time();
214
215 // Output results of move statistics to files
216 mcMoveManager_.output();
218 analyzerManager_.output();
219 }
220
221 // Output results of ramp
222 if (hasRamp()){
223 ramp().output();
224 }
225
226 // Output times for the simulation run
227 Log::file() << std::endl;
228 Log::file() << "nStep " << nStep << std::endl;
229 if (iStep_ != nStep){
230 Log::file() << "nFail Step " << (nStep - iStep_) << std::endl;
231 }
232 Log::file() << "Total run time " << time
233 << " sec" << std::endl;
234 double rStep = double(nStep);
235 Log::file() << "time / nStep " << time / rStep
236 << " sec" << std::endl;
237 Log::file() << "Analyzer run time " << analyzerTime
238 << " sec" << std::endl;
239 Log::file() << std::endl;
240
241 // Output number of times MDE has been solved for the simulation run
242 outputMdeCounter(Log::file());
243
244 // Print McMove acceptance statistics
245 long attempt;
246 long accept;
247 long fail;
248 using namespace std;
249 Log::file() << "Move Statistics:" << endl << endl;
250 Log::file() << setw(20) << left << "Move Name"
251 << setw(10) << right << "Attempted"
252 << setw(10) << right << "Accepted"
253 << setw(13) << right << "AcceptRate"
254 << setw(10) << right << "Failed"
255 << setw(13) << right << "FailRate"
256 << endl;
257 int nMove = mcMoveManager_.size();
258 for (int iMove = 0; iMove < nMove; ++iMove) {
259 attempt = mcMoveManager_[iMove].nAttempt();
260 accept = mcMoveManager_[iMove].nAccept();
261 fail = mcMoveManager_[iMove].nFail();
262 Log::file() << setw(20) << left
263 << mcMoveManager_[iMove].className()
264 << setw(10) << right << attempt
265 << setw(10) << accept
266 << setw(13) << fixed << setprecision(5)
267 << ( attempt == 0 ? 0.0 : double(accept)/double(attempt) )
268 << setw(10) << fail
269 << setw(13) << fixed << setprecision(5)
270 << ( attempt == 0 ? 0.0 : double(fail)/double(attempt) )
271 << endl;
272 }
273 Log::file() << endl;
274
275 }
276
277 /*
278 * Open, read and analyze a trajectory file
279 */
280 template <int D>
281 void McSimulator<D>::analyze(int min, int max,
282 std::string classname,
283 std::string filename)
284 {
285 // Preconditions
286 UTIL_CHECK(min >= 0);
287 UTIL_CHECK(max >= min);
289 UTIL_CHECK(analyzerManager_.size() > 0);
290
291 // Construct TrajectoryReader
292 TrajectoryReader<D>* trajectoryReaderPtr;
293 trajectoryReaderPtr = trajectoryReaderFactory().factory(classname);
294 if (!trajectoryReaderPtr) {
295 std::string message;
296 message = "Invalid TrajectoryReader class name " + classname;
297 UTIL_THROW(message.c_str());
298 }
299
300 // Open trajectory file
301 trajectoryReaderPtr->open(filename);
302 trajectoryReaderPtr->readHeader();
303
304 // Main loop over trajectory frames
305 Timer timer;
306 timer.start();
307 bool hasFrame;
308 hasFrame = trajectoryReaderPtr->readFrame();
309
310 for (iStep_ = 0; iStep_ <= max && hasFrame; ++iStep_) {
311 if (hasFrame) {
312 clearData();
313
314 // Initialize analyzers
315 if (iStep_ == min) {
316 //analyzerManager_.setup();
317 setup(iStep_);
318 }
319
320 // Sample property values only for iStep >= min
321 if (iStep_ >= min) {
322 analyzerManager_.sample(iStep_);
323 }
324 }
325
326 hasFrame = trajectoryReaderPtr->readFrame();
327 }
328 timer.stop();
329 Log::file() << "end main loop" << std::endl;
330 int nFrames = iStep_ - min;
331 trajectoryReaderPtr->close();
332 delete trajectoryReaderPtr;
333
334 // Output results of all analyzers to output files
335 analyzerManager_.output();
336
337 // Output number of frames and times
338 Log::file() << std::endl;
339 Log::file() << "# of frames " << nFrames << std::endl;
340 Log::file() << "run time " << timer.time()
341 << " sec" << std::endl;
342 Log::file() << "time / frame " << timer.time()/double(nFrames)
343 << " sec" << std::endl;
344 Log::file() << std::endl;
345
346 }
347
348 /*
349 * Output McMoveManager timer results.
350 */
351 template<int D>
352 void McSimulator<D>::outputTimers(std::ostream& out)
353 {
354 compressor().outputTimers(out);
355 out << "\n";
356 out << "MC move time contributions:\n";
357 mcMoveManager_.outputTimers(out);
358 }
359
360 /*
361 * Clear all McMoveManager timers.
362 */
363 template<int D>
365 { mcMoveManager_.clearTimers(); }
366
367}
368}
369#endif
Abstract base for periodic output and/or analysis actions.
Monte-Carlo simulation coordinator.
virtual void outputTimers(std::ostream &out)
Output timing results.
McSimulator(System< D > &system)
Constructor.
void setClassName(const char *className)
Set class name string.
virtual void analyze(int min, int max, std::string classname, std::string filename)
Read and analyze a trajectory file.
void simulate(int nStep)
Perform a field theoretic Monte-Carlo simulation.
virtual void clearTimers()
Clear timers.
virtual void readParameters(std::istream &in)
Read parameters for a MC simulation.
Field theoretic simulator (base class).
Definition rpg/System.h:41
System< D > & system()
Get parent system by reference.
Main class for calculations that represent one system.
Definition rpg/System.h:107
Factory for subclasses of TrajectoryReader.
virtual bool readFrame()=0
Read a single frame.
virtual void close()=0
Close the trajectory file.
virtual void open(std::string filename)=0
Open trajectory file and read header, if any.
static std::ostream & file()
Get log ostream by reference.
Definition Log.cpp:57
Wall clock timer.
Definition Timer.h:35
void start(TimePoint begin)
Start timing from an externally supplied time.
Definition Timer.cpp:49
double time()
Return the accumulated time, in seconds.
Definition Timer.cpp:37
void stop(TimePoint end)
Stop the clock at an externally supplied time.
Definition Timer.cpp:73
File containing preprocessor macros for error handling.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
Definition global.h:68
#define UTIL_THROW(msg)
Macro for throwing an Exception, reporting function, file and line number.
Definition global.h:51
Fields, FFTs, and utilities for periodic boundary conditions (CUDA)
Definition CField.cu:12
PSCF package top-level namespace.
Definition param_pc.dox:1
Utility classes for scientific computation.