PSCF v1.2
rpc/fts/compressor/LrCompressor.tpp
1#ifndef RPC_LR_COMPRESSOR_TPP
2#define RPC_LR_COMPRESSOR_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 "LrCompressor.h"
12#include <rpc/System.h>
13#include <rpc/fts/compressor/intra/IntraCorrelation.h>
14#include <util/global.h>
15#include <util/format/Dbl.h>
16#include <pscf/mesh/MeshIterator.h>
17#include <pscf/iterator/NanException.h>
18#include <prdc/crystal/shiftToMinimum.h>
19
20
21namespace Pscf {
22namespace Rpc{
23
24 using namespace Util;
25
26 // Constructor
27 template <int D>
29 : Compressor<D>(system),
30 epsilon_(0.0),
31 itr_(0),
32 maxItr_(0),
33 totalItr_(0),
34 errorType_("rmsResid"),
35 verbose_(0),
36 intra_(system),
37 isIntraCalculated_(false),
38 isAllocated_(false)
39 { setClassName("LrCompressor"); }
40
41 // Destructor
42 template <int D>
45
46 // Read parameters from file
47 template <int D>
48 void LrCompressor<D>::readParameters(std::istream& in)
49 {
50 maxItr_ = 60;
51 read(in, "epsilon", epsilon_);
52 readOptional(in, "maxItr", maxItr_);
53 readOptional(in, "verbose", verbose_);
54 readOptional(in, "errorType", errorType_);
55 }
56
57 // Initialize just before entry to iterative loop.
58 template <int D>
60 {
61 const int nMonomer = system().mixture().nMonomer();
62 IntVec<D> const & dimensions = system().domain().mesh().dimensions();
63 for (int i = 0; i < D; ++i) {
64 if (i < D - 1) {
65 kMeshDimensions_[i] = dimensions[i];
66 } else {
67 kMeshDimensions_[i] = dimensions[i]/2 + 1;
68 }
69 }
70
71 // Allocate memory required by AM algorithm if not done earlier.
72 if (!isAllocated_){
73 resid_.allocate(dimensions);
74 residK_.allocate(dimensions);
75 wFieldTmp_.allocate(nMonomer);
76 intraCorrelationK_.allocate(kMeshDimensions_);
77 for (int i = 0; i < nMonomer; ++i) {
78 wFieldTmp_[i].allocate(dimensions);
79 }
80 isAllocated_ = true;
81 }
82
83 // Compute intraCorrelation
84 if (!isIntraCalculated_){
85 intra_.computeIntraCorrelations(intraCorrelationK_);
86 isIntraCalculated_ = true;
87 }
88 }
89
90 template <int D>
92 {
93 // Initialization and allocate operations on entry to loop.
94 setup();
95 UTIL_CHECK(isAllocated_);
96
97 // Start overall timer
98 timerTotal_.start();
99
100 // Solve MDE
101 timerMDE_.start();
102 system().compute();
103 ++mdeCounter_;
104 timerMDE_.stop();
105
106 // Iterative loop
107 for (itr_ = 0; itr_ < maxItr_; ++itr_) {
108
109 if (verbose_ > 2) {
110 Log::file() << "------------------------------- \n";
111 }
112
113 if (verbose_ > 0){
114 Log::file() << " Iteration " << Int(itr_,5);
115 }
116
117 // Compute residual vector
118 getResidual();
119 double error;
120 try {
121 error = computeError(verbose_);
122 } catch (const NanException&) {
123 Log::file() << ", error = NaN" << std::endl;
124 break; // Exit loop if a NanException is caught
125 }
126 if (verbose_ > 0) {
127 Log::file() << ", error = " << Dbl(error, 15) << std::endl;
128 }
129
130 // Check for convergence
131 if (error < epsilon_) {
132
133 // Successful completion (i.e., converged within tolerance)
134 timerTotal_.stop();
135
136 if (verbose_ > 2) {
137 Log::file() << "-------------------------------\n";
138 }
139 if (verbose_ > 0) {
140 Log::file() << " Converged\n";
141 }
142 if (verbose_ == 2) {
143 Log::file() << "\n";
144 computeError(2);
145 }
146 //mdeCounter_ += itr_;
147 totalItr_ += itr_;
148
149 return 0; // Success
150
151 } else{
152
153 // Not yet converged.
154 updateWFields();
155 timerMDE_.start();
156 system().compute();
157 ++mdeCounter_;
158 timerMDE_.stop();
159
160 }
161
162 }
163
164 // Failure: iteration counter itr reached maxItr without converging
165 timerTotal_.stop();
166
167 Log::file() << "Iterator failed to converge.\n";
168 return 1;
169
170 }
171
172 // Compute the residual for the current system state
173 template <int D>
175 {
176 const int nMonomer = system().mixture().nMonomer();
177 const int meshSize = system().domain().mesh().size();
178
179 // Initialize residuals
180 for (int i = 0 ; i < meshSize; ++i) {
181 resid_[i] = -1.0;
182 }
183
184 // Compute SCF residual vector elements
185 for (int j = 0; j < nMonomer; ++j) {
186 for (int k = 0; k < meshSize; ++k) {
187 resid_[k] += system().c().rgrid(j)[k];
188 }
189 }
190 }
191
192 // update system w field using linear response approximation
193 template <int D>
194 void LrCompressor<D>::updateWFields()
195 {
196 const int nMonomer = system().mixture().nMonomer();
197 const int meshSize = system().domain().mesh().size();
198 const double vMonomer = system().mixture().vMonomer();
199
200 // Convert residual to Fourier Space
201 system().domain().fft().forwardTransform(resid_, residK_);
202 MeshIterator<D> iter;
203 iter.setDimensions(residK_.dftDimensions());
204 for (iter.begin(); !iter.atEnd(); ++iter) {
205 residK_[iter.rank()][0] *= 1.0 / (vMonomer * intraCorrelationK_[iter.rank()]);
206 residK_[iter.rank()][1] *= 1.0 / (vMonomer * intraCorrelationK_[iter.rank()]);
207 }
208
209 // Convert back to real Space (destroys residK_)
210 system().domain().fft().inverseTransformUnsafe(residK_, resid_);
211
212 for (int i = 0; i < nMonomer; i++){
213 for (int k = 0; k < meshSize; k++){
214 wFieldTmp_[i][k] = system().w().rgrid(i)[k] + resid_[k];
215 }
216 }
217 system().setWRGrid(wFieldTmp_);
218 }
219
220 template<int D>
221 void LrCompressor<D>::outputToLog()
222 {}
223
224 template<int D>
225 void LrCompressor<D>::outputTimers(std::ostream& out)
226 {
227 // Output timing results, if requested.
228 double total = timerTotal_.time();
229 out << "\n";
230 out << "LrCompressor time contributions:\n";
231 out << " ";
232 out << "Total" << std::setw(22)<< "Per Iteration"
233 << std::setw(9) << "Fraction" << "\n";
234 out << "MDE solution: "
235 << Dbl(timerMDE_.time(), 9, 3) << " s, "
236 << Dbl(timerMDE_.time()/totalItr_, 9, 3) << " s, "
237 << Dbl(timerMDE_.time()/total, 9, 3) << "\n";
238 out << "\n";
239 }
240
241 template<int D>
243 {
244 timerTotal_.clear();
245 timerMDE_.clear();
246 mdeCounter_ = 0;
247 totalItr_ = 0;
248 }
249
250 template <int D>
251 double LrCompressor<D>::maxAbs(RField<D> const & a)
252 {
253 const int n = a.capacity();
254 double max = 0.0;
255 double value;
256 for (int i = 0; i < n; i++) {
257 value = a[i];
258 if (std::isnan(value)) { // if value is NaN, throw NanException
259 throw NanException("LrCompressor::dotProduct",__FILE__,__LINE__,0);
260 }
261 if (fabs(value) > max) {
262 max = fabs(value);
263 }
264 }
265 return max;
266 }
267
268 // Compute and return inner product of two vectors.
269 template <int D>
270 double LrCompressor<D>::dotProduct(RField<D> const & a,
271 RField<D> const & b)
272 {
273 const int n = a.capacity();
274 UTIL_CHECK(b.capacity() == n);
275 double product = 0.0;
276 for (int i = 0; i < n; i++) {
277 // if either value is NaN, throw NanException
278 if (std::isnan(a[i]) || std::isnan(b[i])) {
279 throw NanException("AmCompressor::dotProduct",__FILE__,__LINE__,0);
280 }
281 product += a[i] * b[i];
282 }
283 return product;
284 }
285
286 // Compute L2 norm of an RField
287 template <int D>
288 double LrCompressor<D>::norm(RField<D> const & a)
289 {
290 double normSq = dotProduct(a, a);
291 return sqrt(normSq);
292 }
293
294 // Compute and return the scalar error
295 template <int D>
296 double LrCompressor<D>::computeError(int verbose)
297 {
298 const int meshSize = system().domain().mesh().size();
299 double error = 0;
300
301 // Find max residual vector element
302 double maxRes = maxAbs(resid_);
303 // Find norm of residual vector
304 double normRes = norm(resid_);
305 // Find root-mean-squared residual element value
306 double rmsRes = normRes/sqrt(meshSize);
307 if (errorType_ == "maxResid") {
308 error = maxRes;
309 } else if (errorType_ == "normResid") {
310 error = normRes;
311 } else if (errorType_ == "rmsResid") {
312 error = rmsRes;
313 } else {
314 UTIL_THROW("Invalid iterator error type in parameter file.");
315 }
316
317 if (verbose > 1) {
318 Log::file() << "\n";
319 Log::file() << "Max Residual = " << Dbl(maxRes,15) << "\n";
320 Log::file() << "Residual Norm = " << Dbl(normRes,15) << "\n";
321 Log::file() << "RMS Residual = " << Dbl(rmsRes,15) << "\n";
322
323 // Check if calculation has diverged (normRes will be NaN)
324 UTIL_CHECK(!std::isnan(normRes));
325 }
326 return error;
327 }
328
329}
330}
331#endif
An IntVec<D, T> is a D-component vector of elements of integer type T.
Definition IntVec.h:27
Exception thrown when not-a-number (NaN) is encountered.
Field of real double precision values on an FFT mesh.
Base class for iterators that impose incompressibility.
void setClassName(const char *className)
Set class name string.
LrCompressor(System< D > &system)
Constructor.
int compress()
Iterate pressure field to obtain partial saddle point.
void outputTimers(std::ostream &out)
Return compressor times contributions.
void readParameters(std::istream &in)
Read all parameters and initialize.
void setup()
Initialize just before entry to iterative loop.
Main class for SCFT or PS-FTS simulation of one system.
Definition rpc/System.h:100
int capacity() const
Return allocated size.
Definition Array.h:159
Wrapper for a double precision number, for formatted ostream output.
Definition Dbl.h:40
Wrapper for an int, for formatted ostream output.
Definition Int.h:37
static std::ostream & file()
Get log ostream by reference.
Definition Log.cpp:57
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
cudaReal maxAbs(DeviceArray< cudaReal > const &in)
Get maximum absolute magnitude of array elements (GPU kernel wrapper).
Definition Reduce.cu:648
cudaReal max(DeviceArray< cudaReal > const &in)
Get maximum of array elements (GPU kernel wrapper).
Definition Reduce.cu:569
PSCF package top-level namespace.
Definition param_pc.dox:1
Utility classes for scientific computation.
float product(float a, float b)
Product for float Data.
Definition product.h:22