PSCF v1.3.3
rpg/fts/compressor/LrCompressor.tpp
1#ifndef RPG_LR_COMPRESSOR_TPP
2#define RPG_LR_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 "LrCompressor.h"
12#include <rpg/solvers/Mixture.h>
13#include <rpg/field/Domain.h>
14#include <rpg/system/System.h>
15#include <prdc/crystal/shiftToMinimum.h>
16#include <prdc/cuda/resources.h>
17#include <pscf/mesh/MeshIterator.h>
18#include <pscf/iterator/NanException.h>
19#include <util/format/Dbl.h>
20#include <util/global.h>
21
22namespace Pscf {
23namespace Rpg{
24
25 using namespace Util;
26
27 /*
28 * Constructor.
29 */
30 template <int D>
32 : Compressor<D>(system),
33 epsilon_(0.0),
34 itr_(0),
35 maxItr_(0),
36 totalItr_(0),
37 verbose_(0),
38 errorType_("rms"),
39 intra_(system),
40 isIntraCalculated_(false),
41 isAllocated_(false)
42 { setClassName("LrCompressor"); }
43
44 /*
45 * Destructor.
46 */
47 template <int D>
50
51 /*
52 * Read parameters from file.
53 */
54 template <int D>
55 void LrCompressor<D>::readParameters(std::istream& in)
56 {
57 // Default values
58 maxItr_ = 100;
59 verbose_ = 0;
60 errorType_ = "rms";
61
62 // Read parameters
63 read(in, "epsilon", epsilon_);
64 readOptional(in, "maxItr", maxItr_);
65 readOptional(in, "verbose", verbose_);
66 readOptional(in, "errorType", errorType_);
67 }
68
69 /*
70 * Initialize just before entry to iterative loop.
71 */
72 template <int D>
74 {
75 const int nMonomer = system().mixture().nMonomer();
76 IntVec<D> const & dimensions = system().domain().mesh().dimensions();
77
78 FFT<D>::computeKMesh(dimensions, kMeshDimensions_, kSize_);
79
80 // Allocate memory required by AM algorithm if not done earlier.
81 if (!isAllocated_){
82 resid_.allocate(dimensions);
83 residK_.allocate(dimensions);
84 wFieldTmp_.allocate(nMonomer);
85 intraCorrelationK_.allocate(kMeshDimensions_);
86 for (int i = 0; i < nMonomer; ++i) {
87 wFieldTmp_[i].allocate(dimensions);
88 }
89 isAllocated_ = true;
90 }
91
92 // Compute intraCorrelation
93 if (!isIntraCalculated_){
94 intra_.computeIntraCorrelations(intraCorrelationK_);
95 isIntraCalculated_ = true;
96 }
97
98 }
99
100 /*
101 * Iteratively adjust pressure field.
102 */
103 template <int D>
105 {
106 // Initialization and allocate operations on entry to loop.
107 setup();
108 UTIL_CHECK(isAllocated_);
109
110 // Start overall timer
111 timerTotal_.start();
112
113 // Solve MDE
114 timerMDE_.start();
115 system().compute();
116 ++mdeCounter_;
117 timerMDE_.stop();
118
119 // Iterative loop
120 for (itr_ = 0; itr_ < maxItr_; ++itr_) {
121
122 if (verbose_ > 2) {
123 Log::file() << "------------------------------- \n";
124 }
125
126 if (verbose_ > 0){
127 Log::file() << " Iteration " << Int(itr_,5);
128 }
129
130 // Compute residual vector
131 getResidual();
132 double error;
133 try {
134 error = computeError(verbose_);
135 } catch (const NanException&) {
136 Log::file() << ", error = NaN" << std::endl;
137 break; // Exit loop if a NanException is caught
138 }
139 if (verbose_ > 0) {
140 Log::file() << ", error = " << Dbl(error, 15) << std::endl;
141 }
142
143 // Check for convergence
144 if (error < epsilon_) {
145
146 // Successful completion (i.e., converged within tolerance)
147 timerTotal_.stop();
148
149 if (verbose_ > 2) {
150 Log::file() << "-------------------------------\n";
151 }
152 if (verbose_ > 0) {
153 Log::file() << " Converged\n";
154 }
155 if (verbose_ == 2) {
156 Log::file() << "\n";
157 computeError(2);
158 }
159 //mdeCounter_ += itr_;
160 totalItr_ += itr_;
161
162 return 0; // Success
163
164 } else{
165
166 // Not yet converged.
167 updateWFields();
168 timerMDE_.start();
169 system().compute();
170 ++mdeCounter_;
171 timerMDE_.stop();
172
173 }
174
175 }
176
177 // Failure: iteration counter itr reached maxItr without converging
178 timerTotal_.stop();
179
180 Log::file() << "Iterator failed to converge.\n";
181 return 1;
182
183 }
184
185 /*
186 * Compute the residual for the current system state.
187 */
188 template <int D>
189 void LrCompressor<D>::getResidual()
190 {
191 const int nMonomer = system().mixture().nMonomer();
192
193 // Initialize resid to c field of species 0 minus 1
194 VecOp::subVS(resid_, system().c().rgrid(0), 1.0);
195
196 // Add other c fields to get SCF residual vector elements
197 for (int i = 1; i < nMonomer; i++) {
198 VecOp::addEqV(resid_, system().c().rgrid(i));
199 }
200 }
201
202 /*
203 * Update system w field using linear response approximation.
204 */
205 template <int D>
206 void LrCompressor<D>::updateWFields(){
207 const int nMonomer = system().mixture().nMonomer();
208 const double vMonomer = system().mixture().vMonomer();
209
210 // Convert residual to Fourier Space
211 system().domain().fft().forwardTransform(resid_, residK_);
212
213 // Compute change in fields using estimated Jacobian
214 VecOp::divEqVc(residK_, intraCorrelationK_, vMonomer);
215
216 // Convert back to real Space (destroys residK_)
217 system().domain().fft().inverseTransformUnsafe(residK_, resid_);
218
219 // Update new fields
220 for (int i = 0; i < nMonomer; i++) {
221 VecOp::addVV(wFieldTmp_[i], system().w().rgrid(i), resid_);
222 }
223
224 // Set system r grid
225 system().w().setRGrid(wFieldTmp_);
226 }
227
228 /*
229 * Output information to a log file (do-nothing implementation).
230 */
231 template<int D>
232 void LrCompressor<D>::outputToLog()
233 {}
234
235 /*
236 * Output timing information to evaluate performance.
237 */
238 template<int D>
239 void LrCompressor<D>::outputTimers(std::ostream& out) const
240 {
241 // Output timing results, if requested.
242 double total = timerTotal_.time();
243 out << "\n";
244 out << "LrCompressor time contributions:\n";
245 out << " ";
246 out << "Total" << std::setw(22)<< "Per Iteration"
247 << std::setw(9) << "Fraction" << "\n";
248 out << "MDE solution: "
249 << Dbl(timerMDE_.time(), 9, 3) << " s, "
250 << Dbl(timerMDE_.time()/totalItr_, 9, 3) << " s, "
251 << Dbl(timerMDE_.time()/total, 9, 3) << "\n";
252 out << "\n";
253 }
254
255 /*
256 * Clear all timers.
257 */
258 template<int D>
260 {
261 timerTotal_.clear();
262 timerMDE_.clear();
263 mdeCounter_ = 0;
264 totalItr_ = 0;
265 }
266
267 template <int D>
268 double LrCompressor<D>::maxAbs(RField<D> const & a)
269 {
270 return Reduce::maxAbs(a);
271 }
272
273 /*
274 * Compute and return inner product of two vectors.
275 */
276 template <int D>
277 double LrCompressor<D>::dotProduct(RField<D> const & a,
278 RField<D> const & b)
279 {
280 UTIL_CHECK(a.capacity() == b.capacity());
281 return Reduce::innerProduct(a, b);
282 }
283
284 /*
285 * Compute the L2 norm of an RField.
286 */
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 /*
295 * Compute and return the scalar error.
296 */
297 template <int D>
298 double LrCompressor<D>::computeError(int verbose)
299 {
300 const int meshSize = system().domain().mesh().size();
301 double error = 0;
302
303 // Find max residual vector element
304 double maxRes = maxAbs(resid_);
305 // Find norm of residual vector
306 double normRes = norm(resid_);
307 // Find root-mean-squared residual element value
308 double rmsRes = normRes/sqrt(meshSize);
309 if (errorType_ == "max") {
310 error = maxRes;
311 } else if (errorType_ == "norm") {
312 error = normRes;
313 } else if (errorType_ == "rms") {
314 error = rmsRes;
315 } else {
316 UTIL_THROW("Invalid iterator error type in parameter file.");
317 }
318
319 if (verbose > 1) {
320 Log::file() << "\n";
321 Log::file() << "Max Residual = " << Dbl(maxRes,15) << "\n";
322 Log::file() << "Residual Norm = " << Dbl(normRes,15) << "\n";
323 Log::file() << "RMS Residual = " << Dbl(rmsRes,15) << "\n";
324
325 // Check if calculation has diverged (normRes will be NaN)
326 UTIL_CHECK(!std::isnan(normRes));
327 }
328 return error;
329 }
330
331}
332}
333#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.
static void computeKMesh(IntVec< D > const &rMeshDimensions, IntVec< D > &kMeshDimensions, int &kSize)
Compute dimensions and size of k-space mesh for DFT of real data.
Definition cpu/FFT.tpp:262
Field of real double precision values on an FFT mesh.
Definition cpu/RField.h:29
int mdeCounter_
Count how many times MDE has been solved.
Compressor()
Default constructor.
System< D > const & system() const
Return const reference to parent system.
void setup()
Initialize just before entry to iterative loop.
void readParameters(std::istream &in)
Read all parameters and initialize.
ScalarParam< Type > & read(std::istream &in, const char *label, Type &value)
Add and read a new required ScalarParam < Type > object.
void setClassName(const char *className)
Set class name string.
ScalarParam< Type > & readOptional(std::istream &in, const char *label, Type &value)
Add and read a new optional ScalarParam < Type > object.
int compress()
Iterate pressure field to obtain partial saddle point.
LrCompressor(System< D > &system)
Constructor.
void outputTimers(std::ostream &out) const
Return compressor times contributions.
Main class, representing a complete physical system.
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:59
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:49
double maxAbs(Array< double > const &in)
Get maximum absolute magnitude of array elements .
Definition Reduce.cpp:34
double innerProduct(Array< double > const &a, Array< double > const &b)
Compute inner product of two real arrays .
Definition Reduce.cpp:94
void addEqV(Array< double > &a, Array< double > const &b)
Vector addition in-place, a[i] += b[i].
Definition VecOp.cpp:199
void addVV(Array< double > &a, Array< double > const &b, Array< double > const &c)
Vector addition, a[i] = b[i] + c[i].
Definition VecOp.cpp:50
void subVS(Array< double > &a, Array< double > const &b, double c)
Vector subtraction, a[i] = b[i] - c.
Definition VecOp.cpp:96
void divEqVc(DeviceArray< cudaComplex > &a, DeviceArray< cudaReal > const &b, cudaReal const c)
Vector division in-place w/ coeff., a[i] /= (b[i] * c), kernel wrapper.
Definition VecOpMisc.cu:377
SCFT and PS-FTS with real periodic fields (GPU)
PSCF package top-level namespace.