PSCF v1.3.2
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/system/System.h>
13#include <prdc/crystal/shiftToMinimum.h>
14#include <prdc/cuda/resources.h>
15#include <pscf/mesh/MeshIterator.h>
16#include <pscf/iterator/NanException.h>
17#include <util/format/Dbl.h>
18#include <util/global.h>
19
20namespace Pscf {
21namespace Rpg{
22
23 using namespace Util;
24
25 /*
26 * Constructor.
27 */
28 template <int D>
30 : Compressor<D>(system),
31 epsilon_(0.0),
32 itr_(0),
33 maxItr_(0),
34 totalItr_(0),
35 verbose_(0),
36 errorType_("rms"),
37 intra_(system),
38 isIntraCalculated_(false),
39 isAllocated_(false)
40 { setClassName("LrCompressor"); }
41
42 /*
43 * Destructor.
44 */
45 template <int D>
48
49 /*
50 * Read parameters from file.
51 */
52 template <int D>
53 void LrCompressor<D>::readParameters(std::istream& in)
54 {
55 // Default values
56 maxItr_ = 100;
57 verbose_ = 0;
58 errorType_ = "rms";
59
60 // Read parameters
61 read(in, "epsilon", epsilon_);
62 readOptional(in, "maxItr", maxItr_);
63 readOptional(in, "verbose", verbose_);
64 readOptional(in, "errorType", errorType_);
65 }
66
67 /*
68 * Initialize just before entry to iterative loop.
69 */
70 template <int D>
72 {
73 const int nMonomer = system().mixture().nMonomer();
74 IntVec<D> const & dimensions = system().domain().mesh().dimensions();
75
76 FFT<D>::computeKMesh(dimensions, kMeshDimensions_, kSize_);
77
78 // Allocate memory required by AM algorithm if not done earlier.
79 if (!isAllocated_){
80 resid_.allocate(dimensions);
81 residK_.allocate(dimensions);
82 wFieldTmp_.allocate(nMonomer);
83 intraCorrelationK_.allocate(kMeshDimensions_);
84 for (int i = 0; i < nMonomer; ++i) {
85 wFieldTmp_[i].allocate(dimensions);
86 }
87 isAllocated_ = true;
88 }
89
90 // Compute intraCorrelation
91 if (!isIntraCalculated_){
92 intra_.computeIntraCorrelations(intraCorrelationK_);
93 isIntraCalculated_ = true;
94 }
95
96 }
97
98 /*
99 * Iteratively adjust pressure field.
100 */
101 template <int D>
103 {
104 // Initialization and allocate operations on entry to loop.
105 setup();
106 UTIL_CHECK(isAllocated_);
107
108 // Start overall timer
109 timerTotal_.start();
110
111 // Solve MDE
112 timerMDE_.start();
113 system().compute();
114 ++mdeCounter_;
115 timerMDE_.stop();
116
117 // Iterative loop
118 for (itr_ = 0; itr_ < maxItr_; ++itr_) {
119
120 if (verbose_ > 2) {
121 Log::file() << "------------------------------- \n";
122 }
123
124 if (verbose_ > 0){
125 Log::file() << " Iteration " << Int(itr_,5);
126 }
127
128 // Compute residual vector
129 getResidual();
130 double error;
131 try {
132 error = computeError(verbose_);
133 } catch (const NanException&) {
134 Log::file() << ", error = NaN" << std::endl;
135 break; // Exit loop if a NanException is caught
136 }
137 if (verbose_ > 0) {
138 Log::file() << ", error = " << Dbl(error, 15) << std::endl;
139 }
140
141 // Check for convergence
142 if (error < epsilon_) {
143
144 // Successful completion (i.e., converged within tolerance)
145 timerTotal_.stop();
146
147 if (verbose_ > 2) {
148 Log::file() << "-------------------------------\n";
149 }
150 if (verbose_ > 0) {
151 Log::file() << " Converged\n";
152 }
153 if (verbose_ == 2) {
154 Log::file() << "\n";
155 computeError(2);
156 }
157 //mdeCounter_ += itr_;
158 totalItr_ += itr_;
159
160 return 0; // Success
161
162 } else{
163
164 // Not yet converged.
165 updateWFields();
166 timerMDE_.start();
167 system().compute();
168 ++mdeCounter_;
169 timerMDE_.stop();
170
171 }
172
173 }
174
175 // Failure: iteration counter itr reached maxItr without converging
176 timerTotal_.stop();
177
178 Log::file() << "Iterator failed to converge.\n";
179 return 1;
180
181 }
182
183 /*
184 * Compute the residual for the current system state.
185 */
186 template <int D>
187 void LrCompressor<D>::getResidual()
188 {
189 const int nMonomer = system().mixture().nMonomer();
190
191 // Initialize resid to c field of species 0 minus 1
192 VecOp::subVS(resid_, system().c().rgrid(0), 1.0);
193
194 // Add other c fields to get SCF residual vector elements
195 for (int i = 1; i < nMonomer; i++) {
196 VecOp::addEqV(resid_, system().c().rgrid(i));
197 }
198 }
199
200 /*
201 * Update system w field using linear response approximation.
202 */
203 template <int D>
204 void LrCompressor<D>::updateWFields(){
205 const int nMonomer = system().mixture().nMonomer();
206 const double vMonomer = system().mixture().vMonomer();
207
208 // Convert residual to Fourier Space
209 system().domain().fft().forwardTransform(resid_, residK_);
210
211 // Compute change in fields using estimated Jacobian
212 VecOp::divEqVc(residK_, intraCorrelationK_, vMonomer);
213
214 // Convert back to real Space (destroys residK_)
215 system().domain().fft().inverseTransformUnsafe(residK_, resid_);
216
217 // Update new fields
218 for (int i = 0; i < nMonomer; i++) {
219 VecOp::addVV(wFieldTmp_[i], system().w().rgrid(i), resid_);
220 }
221
222 // Set system r grid
223 system().w().setRGrid(wFieldTmp_);
224 }
225
226 /*
227 * Output information to a log file (do-nothing implementation).
228 */
229 template<int D>
230 void LrCompressor<D>::outputToLog()
231 {}
232
233 /*
234 * Output timing information to evaluate performance.
235 */
236 template<int D>
237 void LrCompressor<D>::outputTimers(std::ostream& out) const
238 {
239 // Output timing results, if requested.
240 double total = timerTotal_.time();
241 out << "\n";
242 out << "LrCompressor time contributions:\n";
243 out << " ";
244 out << "Total" << std::setw(22)<< "Per Iteration"
245 << std::setw(9) << "Fraction" << "\n";
246 out << "MDE solution: "
247 << Dbl(timerMDE_.time(), 9, 3) << " s, "
248 << Dbl(timerMDE_.time()/totalItr_, 9, 3) << " s, "
249 << Dbl(timerMDE_.time()/total, 9, 3) << "\n";
250 out << "\n";
251 }
252
253 /*
254 * Clear all timers.
255 */
256 template<int D>
258 {
259 timerTotal_.clear();
260 timerMDE_.clear();
261 mdeCounter_ = 0;
262 totalItr_ = 0;
263 }
264
265 template <int D>
266 double LrCompressor<D>::maxAbs(RField<D> const & a)
267 {
268 return Reduce::maxAbs(a);
269 }
270
271 /*
272 * Compute and return inner product of two vectors.
273 */
274 template <int D>
275 double LrCompressor<D>::dotProduct(RField<D> const & a,
276 RField<D> const & b)
277 {
278 UTIL_CHECK(a.capacity() == b.capacity());
279 return Reduce::innerProduct(a, b);
280 }
281
282 /*
283 * Compute the L2 norm of an RField.
284 */
285 template <int D>
286 double LrCompressor<D>::norm(RField<D> const & a)
287 {
288 double normSq = dotProduct(a, a);
289 return sqrt(normSq);
290 }
291
292 /*
293 * Compute and return the scalar error.
294 */
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_ == "max") {
308 error = maxRes;
309 } else if (errorType_ == "norm") {
310 error = normRes;
311 } else if (errorType_ == "rms") {
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.
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.