PSCF v1.2
rpc/fts/compressor/LrAmCompressor.tpp
1#ifndef RPC_LR_POST_AM_COMPRESSOR_TPP
2#define RPC_LR_POST_AM_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 "LrAmCompressor.h"
12#include <rpc/System.h>
13#include <rpc/fts/compressor/intra/IntraCorrelation.h>
14#include <pscf/mesh/MeshIterator.h>
15#include <util/global.h>
16
17namespace Pscf {
18namespace Rpc{
19
20 using namespace Util;
21
22 // Constructor
23 template <int D>
25 : Compressor<D>(system),
26 intra_(system),
27 isIntraCalculated_(false),
28 isAllocated_(false)
29 { setClassName("LrAmCompressor"); }
30
31 // Destructor
32 template <int D>
35
36 // Read parameters from file
37 template <int D>
38 void LrAmCompressor<D>::readParameters(std::istream& in)
39 {
40 // Call parent class readParameters
41 AmIteratorTmpl<Compressor<D>, DArray<double> >::readParameters(in);
42 AmIteratorTmpl<Compressor<D>, DArray<double> >::readErrorType(in);
43
44 }
45
46 // Initialize just before entry to iterative loop.
47 template <int D>
48 void LrAmCompressor<D>::setup(bool isContinuation)
49 {
50
51 // Allocate memory required by AM algorithm if not done earlier.
52 AmIteratorTmpl<Compressor<D>, DArray<double> >::setup(isContinuation);
53
54 const int nMonomer = system().mixture().nMonomer();
55 const int meshSize = system().domain().mesh().size();
56 IntVec<D> const & dimensions = system().domain().mesh().dimensions();
57
58 // Compute Fourier space kMeshDimensions_
59 for (int i = 0; i < D; ++i) {
60 if (i < D - 1) {
61 kMeshDimensions_[i] = dimensions[i];
62 } else {
63 kMeshDimensions_[i] = dimensions[i]/2 + 1;
64 }
65 }
66
67 // Allocate memory required by compressor, if not done previously
68 if (!isAllocated_){
69 newBasis_.allocate(meshSize);
70 w0_.allocate(nMonomer);
71 wFieldTmp_.allocate(nMonomer);
72 resid_.allocate(dimensions);
73 residK_.allocate(dimensions);
74 intraCorrelationK_.allocate(kMeshDimensions_);
75 for (int i = 0; i < nMonomer; ++i) {
76 w0_[i].allocate(dimensions);
77 wFieldTmp_[i].allocate(dimensions);
78 }
79
80 isAllocated_ = true;
81 }
82
83 // Compute intraCorrelation
84 if (!isIntraCalculated_){
85 intra_.computeIntraCorrelations(intraCorrelationK_);
86 isIntraCalculated_ = true;
87 }
88
89 // Store initial values of monomer chemical potential fields in w0_
90 for (int i = 0; i < nMonomer; ++i) {
91 for (int j = 0; j< meshSize; ++j){
92 w0_[i][j] = system().w().rgrid(i)[j];
93 }
94 }
95
96 }
97
98 /*
99 * Apply the Anderson-Mixing algorithm (main function).
100 */
101 template <int D>
103 {
104 int solve = AmIteratorTmpl<Compressor<D>, DArray<double> >::solve();
105 return solve;
106 }
107
108 /*
109 * Assign one array to another.
110 */
111 template <int D>
113 DArray<double> const & b)
114 { a = b; }
115
116 /*
117 * Compute and return inner product of two vectors.
118 */
119 template <int D>
120 double
121 LrAmCompressor<D>::dotProduct(DArray<double> const & a,
122 DArray<double> const & b)
123 {
124 const int n = a.capacity();
125 UTIL_CHECK(b.capacity() == n);
126 double product = 0.0;
127 for (int i = 0; i < n; i++) {
128 // if either value is NaN, throw NanException
129 if (std::isnan(a[i]) || std::isnan(b[i])) {
130 throw NanException("LrAmCompressor::dotProduct",
131 __FILE__,__LINE__,0);
132 }
133 product += a[i] * b[i];
134 }
135 return product;
136 }
137
138 /*
139 * Compute and return the maximum element of a vector.
140 */
141 template <int D>
142 double LrAmCompressor<D>::maxAbs(DArray<double> const & a)
143 {
144 const int n = a.capacity();
145 double max = 0.0;
146 double value;
147 for (int i = 0; i < n; i++) {
148 value = a[i];
149 if (std::isnan(value)) { // if value is NaN, throw NanException
150 throw NanException("LrAmCompressor::dotProduct",
151 __FILE__, __LINE__, 0);
152 }
153 if (fabs(value) > max) {
154 max = fabs(value);
155 }
156 }
157 return max;
158 }
159
160 /*
161 * Update basis by adding one new basis vector.
162 */
163 template <int D>
164 void
165 LrAmCompressor<D>::updateBasis(
166 RingBuffer< DArray<double> > & basis,
167 RingBuffer< DArray<double> > const & hists)
168 {
169 // Make sure at least two histories are stored
170 UTIL_CHECK(hists.size() >= 2);
171
172 const int n = hists[0].capacity();
173
174 // New basis vector is difference between two most recent states
175 for (int i = 0; i < n; i++) {
176 newBasis_[i] = hists[0][i] - hists[1][i];
177 }
178
179 basis.append(newBasis_);
180 }
181
182 template <int D>
183 void
184 LrAmCompressor<D>::addHistories(DArray<double>& trial,
185 RingBuffer<DArray<double> > const & basis,
186 DArray<double> coeffs,
187 int nHist)
188 {
189 int n = trial.capacity();
190 for (int i = 0; i < nHist; i++) {
191 for (int j = 0; j < n; j++) {
192 // Not clear on the origin of the -1 factor
193 trial[j] += coeffs[i] * -1 * basis[i][j];
194 }
195 }
196 }
197
198 /*
199 * Returns mixing value parameter lambda = 1.0
200 */
201 template<int D>
202 double LrAmCompressor<D>::computeLambda(double r)
203 { return 1.0; }
204
205 /*
206 * Correction step (second step of Anderson mixing)
207 *
208 * This LrAM algorithm uses a quasi-Newton correction step with an
209 * approximate Jacobian given by the Jacobian in a homogeneous state.
210 */
211 template <int D>
212 void
213 LrAmCompressor<D>::addPredictedError(DArray<double>& fieldTrial,
214 DArray<double> const & resTrial,
215 double lambda)
216 {
217 // Local constants
218 const int n = fieldTrial.capacity();
219 const double vMonomer = system().mixture().vMonomer();
220
221 // Copy DArray<double> resTrial into RField<D> resid_
222 // Allows use of FFT functions that take an RField container
223 for (int i = 0 ; i < n; ++i) {
224 resid_[i] = resTrial[i];
225 }
226
227 // Fourier transform r-grid residual resid_ to k-grid form residK_
228 system().domain().fft().forwardTransform(resid_, residK_);
229
230 // Compute update on a k-grid, using quasi-Newton algorithm
231 double factor;
232 MeshIterator<D> iter;
233 iter.setDimensions(kMeshDimensions_);
234 for (iter.begin(); !iter.atEnd(); ++iter) {
235 factor = 1.0 / (vMonomer * intraCorrelationK_[iter.rank()]);
236 residK_[iter.rank()][0] *= factor;
237 residK_[iter.rank()][1] *= factor;
238 }
239 // Field residK_ now stores k-grid update rather than residual
240
241 // Convert field update back to r-grid (real space) form
242 // On return, resid_ contains r-grid field update, residK_ is destroyed
243 system().domain().fft().inverseTransformUnsafe(residK_, resid_);
244
245 // Add update to obtain new fieldTrial
246 for (int i = 0; i < n; i++) {
247 fieldTrial[i] += lambda * resid_[i];
248 }
249 }
250
251 /*
252 * Does the associated system have initialized w fields?
253 */
254 template <int D>
255 bool LrAmCompressor<D>::hasInitialGuess()
256 { return system().w().hasData(); }
257
258 /*
259 * Compute and return the number of elements in a field vector.
260 */
261 template <int D>
262 int LrAmCompressor<D>::nElements()
263 { return system().domain().mesh().size(); }
264
265 /*
266 * Get the current field variable from the system.
267 *
268 * The field variable is the change in the Lagrange multiplier field.
269 * relative to that used in initial array of monomer fields, w0_.
270 */
271 template <int D>
272 void LrAmCompressor<D>::getCurrent(DArray<double>& curr)
273 {
274 // Straighten out fields into linear arrays
275 const int meshSize = system().domain().mesh().size();
276 const DArray< RField<D> > * currSys = &system().w().rgrid();
277
278 /*
279 * The field that we are adjusting is the Langrange multiplier field
280 * with number of grid pts components. The current value is the
281 * difference between w and w0_ for the first monomer (any monomer
282 * should give the same answer)
283 */
284 for (int i = 0; i < meshSize; i++){
285 curr[i] = (*currSys)[0][i] - w0_[0][i];
286 }
287
288 }
289
290 /*
291 * Perform the main system computation (solve the MDE).
292 */
293 template <int D>
294 void LrAmCompressor<D>::evaluate()
295 {
296 system().compute();
297 ++mdeCounter_;
298 }
299
300 /*
301 * Compute the residual vector for the current system state
302 */
303 template <int D>
304 void LrAmCompressor<D>::getResidual(DArray<double>& resid)
305 {
306 const int n = nElements();
307 const int nMonomer = system().mixture().nMonomer();
308 const int meshSize = system().domain().mesh().size();
309
310 // Initialize residuals
311 for (int i = 0 ; i < n; ++i) {
312 resid[i] = -1.0;
313 }
314
315 // Compute SCF residual vector elements
316 for (int j = 0; j < nMonomer; ++j) {
317 for (int k = 0; k < meshSize; ++k) {
318 resid[k] += system().c().rgrid(j)[k];
319 }
320 }
321
322 }
323
324 /*
325 * Update the w field values stored in the system
326 */
327 template <int D>
328 void LrAmCompressor<D>::update(DArray<double>& newGuess)
329 {
330 // Local constant values
331 const int nMonomer = system().mixture().nMonomer();
332 const int meshSize = system().domain().mesh().size();
333
334 // Locally update monomer w fields in array wFieldTmp_
335 // Set wFieldTmp = w0_ + newGuess
336 for (int i = 0; i < nMonomer; i++){
337 for (int k = 0; k < meshSize; k++){
338 wFieldTmp_[i][k] = w0_[i][k] + newGuess[k];
339 }
340 }
341
342 // Update w fields in the system WContainer
343 system().setWRGrid(wFieldTmp_);
344 }
345
346 template<int D>
347 void LrAmCompressor<D>::outputToLog()
348 {}
349
350 /*
351 * Output timing results to an output stream (file)
352 */
353 template<int D>
354 void LrAmCompressor<D>::outputTimers(std::ostream& out)
355 {
356 // Output timing results, if requested.
357 out << "\n";
358 out << "LrAmCompressor time contributions:\n";
359 AmIteratorTmpl<Compressor<D>, DArray<double> >::outputTimers(out);
360 }
361
362 /*
363 * Clear all timers and the MDE counter
364 */
365 template<int D>
367 {
369 mdeCounter_ = 0;
370 }
371
372}
373}
374#endif
Template for Anderson mixing iterator algorithm.
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.
Base class for iterators that impose incompressibility.
Anderson Mixing compressor with linear-response mixing step.
void readParameters(std::istream &in)
Read all parameters and initialize.
LrAmCompressor(System< D > &system)
Constructor.
void outputTimers(std::ostream &out)
Return compressor times contributions.
void clearTimers()
Clear all timers (reset accumulated time to zero).
void setup(bool isContinuation)
Initialize just before entry to iterative loop.
int compress()
Compress to obtain partial saddle point w+.
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
Dynamically allocatable contiguous array template.
Class for storing history of previous values in an array.
Definition RingBuffer.h:27
File containing preprocessor macros for error handling.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
Definition global.h:68
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