PSCF v1.3.2
rpc/fts/compressor/LrAmCompressor.tpp
1#ifndef RPC_LR_AM_COMPRESSOR_TPP
2#define RPC_LR_AM_COMPRESSOR_TPP
3
4/*
5* PSCF - Polymer Self-Consistent Field
6*
7* Copyright 2015 - 2025, 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/System.h>
13#include <pscf/mesh/MeshIterator.h>
14#include <util/global.h>
15
16namespace Pscf {
17namespace Rpc{
18
19 using namespace Util;
20
21 // Public member functions
22
23 /*
24 * Constructor.
25 */
26 template <int D>
28 : Compressor<D>(system),
29 intra_(system),
30 isIntraCalculated_(false),
31 isAllocated_(false)
32 { ParamComposite::setClassName("LrAmCompressor"); }
33
34 /*
35 * Destructor.
36 */
37 template <int D>
40
41 /*
42 * Read parameters from file.
43 */
44 template <int D>
45 void LrAmCompressor<D>::readParameters(std::istream& in)
46 {
47 // Default values
48 AmTmpl::maxItr_ = 100;
50 AmTmpl::errorType_ = "rms";
51
52 // Call parent read methods
55
56 }
57
58 /*
59 * Initialize just before entry to iterative loop.
60 */
61 template <int D>
62 void LrAmCompressor<D>::setup(bool isContinuation)
63 {
64
65 // Allocate memory required by AM algorithm if not done earlier.
67
68 const int nMonomer = system().mixture().nMonomer();
69 const int meshSize = system().domain().mesh().size();
70 IntVec<D> const & dimensions = system().domain().mesh().dimensions();
71
72 // Compute Fourier space kMeshDimensions_
73 for (int i = 0; i < D; ++i) {
74 if (i < D - 1) {
75 kMeshDimensions_[i] = dimensions[i];
76 } else {
77 kMeshDimensions_[i] = dimensions[i]/2 + 1;
78 }
79 }
80
81 // Allocate memory required by compressor, if not done previously
82 if (!isAllocated_){
83 newBasis_.allocate(meshSize);
84 w0_.allocate(nMonomer);
85 wFieldTmp_.allocate(nMonomer);
86 resid_.allocate(dimensions);
87 residK_.allocate(dimensions);
88 intraCorrelationK_.allocate(kMeshDimensions_);
89 for (int i = 0; i < nMonomer; ++i) {
90 w0_[i].allocate(dimensions);
91 wFieldTmp_[i].allocate(dimensions);
92 }
93
94 isAllocated_ = true;
95 }
96
97 // Compute intraCorrelation
98 if (!isIntraCalculated_){
99 intra_.computeIntraCorrelations(intraCorrelationK_);
100 isIntraCalculated_ = true;
101 }
102
103 // Store initial values of monomer chemical potential fields in w0_
104 for (int i = 0; i < nMonomer; ++i) {
105 for (int j = 0; j< meshSize; ++j){
106 w0_[i][j] = system().w().rgrid(i)[j];
107 }
108 }
109
110 }
111
112 /*
113 * Apply the Anderson-Mixing algorithm (main function).
114 */
115 template <int D>
121
122 /*
123 * Output timing results to an output stream (file)
124 */
125 template<int D>
126 void LrAmCompressor<D>::outputTimers(std::ostream& out) const
127 {
128 out << "\n";
129 out << "LrAmCompressor time contributions:\n";
131 }
132
133 /*
134 * Clear all timers and the MDE counter
135 */
136 template<int D>
142
143 // Private virtual function that interact with parent system
144
145 /*
146 * Correction step (second step of Anderson mixing)
147 *
148 * This LrAM algorithm uses a quasi-Newton correction step with an
149 * approximate Jacobian given by the Jacobian in a homogeneous state.
150 */
151 template <int D>
152 void
153 LrAmCompressor<D>::addCorrection(DArray<double>& fieldTrial,
154 DArray<double> const & resTrial)
155 {
156 // Local constants
157 const int n = fieldTrial.capacity();
158 const double vMonomer = system().mixture().vMonomer();
159
160 // Copy DArray<double> resTrial into RField<D> resid_
161 // Allows use of FFT functions that take an RField container
162 for (int i = 0 ; i < n; ++i) {
163 resid_[i] = resTrial[i];
164 }
165
166 // Fourier transform r-grid residual resid_ to k-grid form residK_
167 system().domain().fft().forwardTransform(resid_, residK_);
168
169 // Compute update on a k-grid, using quasi-Newton algorithm
170 double factor;
171 MeshIterator<D> iter;
172 iter.setDimensions(kMeshDimensions_);
173 for (iter.begin(); !iter.atEnd(); ++iter) {
174 factor = 1.0 / (vMonomer * intraCorrelationK_[iter.rank()]);
175 residK_[iter.rank()][0] *= factor;
176 residK_[iter.rank()][1] *= factor;
177 }
178 // Field residK_ now stores k-grid update rather than residual
179
180 // Convert field update back to r-grid (real space) form
181 // On return, resid_ contains r-grid field update, residK_ is destroyed
182 system().domain().fft().inverseTransformUnsafe(residK_, resid_);
183
184 // Add update to obtain new fieldTrial
185 for (int i = 0; i < n; i++) {
186 fieldTrial[i] += resid_[i];
187 }
188 }
189
190 // Private virtual functions that interact with parent System
191
192 /*
193 * Compute and return the number of elements in a field vector.
194 */
195 template <int D>
196 int LrAmCompressor<D>::nElements()
197 { return system().domain().mesh().size(); }
198
199 /*
200 * Does the associated system have initialized w fields?
201 */
202 template <int D>
203 bool LrAmCompressor<D>::hasInitialGuess()
204 { return system().w().hasData(); }
205
206 /*
207 * Get the current field variable from the system.
208 *
209 * The field variable is the change in the Lagrange multiplier field.
210 * relative to that used in initial array of monomer fields, w0_.
211 */
212 template <int D>
213 void LrAmCompressor<D>::getCurrent(DArray<double>& curr)
214 {
215 // Straighten out fields into linear arrays
216 const int meshSize = system().domain().mesh().size();
217 const DArray< RField<D> > * currSys = &system().w().rgrid();
218
219 /*
220 * The field that we are adjusting is the Langrange multiplier field
221 * with number of grid pts components. The current value is the
222 * difference between w and w0_ for the first monomer (any monomer
223 * should give the same answer)
224 */
225 for (int i = 0; i < meshSize; i++){
226 curr[i] = (*currSys)[0][i] - w0_[0][i];
227 }
228
229 }
230
231 /*
232 * Perform the main system computation (solve the MDE).
233 */
234 template <int D>
235 void LrAmCompressor<D>::evaluate()
236 {
237 system().compute();
239 }
240
241 /*
242 * Compute the residual vector for the current system state
243 */
244 template <int D>
245 void LrAmCompressor<D>::getResidual(DArray<double>& resid)
246 {
247 const int n = nElements();
248 const int nMonomer = system().mixture().nMonomer();
249 const int meshSize = system().domain().mesh().size();
250
251 // Initialize residuals
252 for (int i = 0 ; i < n; ++i) {
253 resid[i] = -1.0;
254 }
255
256 // Compute SCF residual vector elements
257 for (int j = 0; j < nMonomer; ++j) {
258 for (int k = 0; k < meshSize; ++k) {
259 resid[k] += system().c().rgrid(j)[k];
260 }
261 }
262
263 }
264
265 /*
266 * Update the w field values stored in the system
267 */
268 template <int D>
269 void LrAmCompressor<D>::update(DArray<double>& newGuess)
270 {
271 // Local constant values
272 const int nMonomer = system().mixture().nMonomer();
273 const int meshSize = system().domain().mesh().size();
274
275 // Locally update monomer w fields in array wFieldTmp_
276 // Set wFieldTmp = w0_ + newGuess
277 for (int i = 0; i < nMonomer; i++){
278 for (int k = 0; k < meshSize; k++){
279 wFieldTmp_[i][k] = w0_[i][k] + newGuess[k];
280 }
281 }
282
283 // Update w fields in the system WContainer
284 system().w().setRGrid(wFieldTmp_);
285 }
286
287 /*
288 * Output results to log file (do-nothing implementation).
289 */
290 template<int D>
291 void LrAmCompressor<D>::outputToLog()
292 {}
293
294}
295}
296#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
Iterator over points in a Mesh<D>.
int rank() const
Get the rank of current element.
void begin()
Set iterator to the first point in the mesh.
bool atEnd() const
Is this the end (i.e., one past the last point)?
void setDimensions(const IntVec< D > &dimensions)
Set the grid dimensions in all directions.
Base class for iterators that impose incompressibility.
int mdeCounter_
Count how many times MDE has been solved.
void readParameters(std::istream &in) override
Read all parameters and initialize.
LrAmCompressor(System< D > &system)
Constructor.
int compress() override
Compress to obtain partial saddle point w+.
void outputTimers(std::ostream &out) const override
Return compressor times contributions.
void setup(bool isContinuation) override
Initialize just before entry to iterative loop.
void clearTimers() override
Clear all timers (reset accumulated time to zero).
Main class, representing a complete physical system.
int capacity() const
Return allocated size.
Definition Array.h:159
Dynamically allocatable contiguous array template.
Definition DArray.h:32
void setClassName(const char *className)
Set class name string.
File containing preprocessor macros for error handling.
Real periodic fields, SCFT and PS-FTS (CPU).
Definition param_pc.dox:2
PSCF package top-level namespace.