PSCF v1.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 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 <rpg/System.h>
13#include <rpg/fts/compressor/intra/IntraCorrelation.h>
14#include <prdc/crystal/shiftToMinimum.h>
15#include <prdc/cuda/resources.h>
16#include <pscf/mesh/MeshIterator.h>
17#include <pscf/iterator/NanException.h>
18#include <util/global.h>
19#include <util/format/Dbl.h>
20
21
22namespace Pscf {
23namespace Rpg{
24
25 using namespace Util;
26
27 // Constructor
28 template <int D>
30 : Compressor<D>(system),
31 epsilon_(0.0),
32 itr_(0),
33 maxItr_(0),
34 totalItr_(0),
35 errorType_("rmsResid"),
36 verbose_(0),
37 intra_(system),
38 isIntraCalculated_(false),
39 isAllocated_(false)
40 { setClassName("LrCompressor"); }
41
42 // Destructor
43 template <int D>
46
47 // Read parameters from file
48 template <int D>
49 void LrCompressor<D>::readParameters(std::istream& in)
50 {
51 // Call parent class readParameters
52 maxItr_ = 60;
53 read(in, "epsilon", epsilon_);
54 readOptional(in, "maxItr", maxItr_);
55 readOptional(in, "verbose", verbose_);
56 readOptional(in, "errorType", errorType_);
57 }
58
59 // Initialize just before entry to iterative loop.
60 template <int D>
62 {
63 const int nMonomer = system().mixture().nMonomer();
64 IntVec<D> const & dimensions = system().mesh().dimensions();
65 for (int i = 0; i < D; ++i) {
66 if (i < D - 1) {
67 kMeshDimensions_[i] = dimensions[i];
68 } else {
69 kMeshDimensions_[i] = dimensions[i]/2 + 1;
70 }
71 }
72
73 // Compute number of points in k-space grid
74 kSize_ = 1;
75 for (int i = 0; i < D; ++i) {
76 kSize_ *= kMeshDimensions_[i];
77 }
78
79 // Allocate memory required by AM algorithm if not done earlier.
80 if (!isAllocated_){
81 resid_.allocate(dimensions);
82 residK_.allocate(dimensions);
83 wFieldTmp_.allocate(nMonomer);
84 intraCorrelationK_.allocate(kMeshDimensions_);
85 for (int i = 0; i < nMonomer; ++i) {
86 wFieldTmp_[i].allocate(dimensions);
87 }
88 isAllocated_ = true;
89 }
90
91 // Compute intraCorrelation
92 if (!isIntraCalculated_){
93 intra_.computeIntraCorrelations(intraCorrelationK_);
94 isIntraCalculated_ = true;
95 }
96
97 }
98
99 template <int D>
101 {
102 // Initialization and allocate operations on entry to loop.
103 setup();
104 UTIL_CHECK(isAllocated_);
105
106 // Start overall timer
107 timerTotal_.start();
108
109 // Solve MDE
110 timerMDE_.start();
111 system().compute();
112 ++mdeCounter_;
113 timerMDE_.stop();
114
115 // Iterative loop
116 for (itr_ = 0; itr_ < maxItr_; ++itr_) {
117
118 if (verbose_ > 2) {
119 Log::file() << "------------------------------- \n";
120 }
121
122 if (verbose_ > 0){
123 Log::file() << " Iteration " << Int(itr_,5);
124 }
125
126 // Compute residual vector
127 getResidual();
128 double error;
129 try {
130 error = computeError(verbose_);
131 } catch (const NanException&) {
132 Log::file() << ", error = NaN" << std::endl;
133 break; // Exit loop if a NanException is caught
134 }
135 if (verbose_ > 0) {
136 Log::file() << ", error = " << Dbl(error, 15) << std::endl;
137 }
138
139 // Check for convergence
140 if (error < epsilon_) {
141
142 // Successful completion (i.e., converged within tolerance)
143 timerTotal_.stop();
144
145 if (verbose_ > 2) {
146 Log::file() << "-------------------------------\n";
147 }
148 if (verbose_ > 0) {
149 Log::file() << " Converged\n";
150 }
151 if (verbose_ == 2) {
152 Log::file() << "\n";
153 computeError(2);
154 }
155 //mdeCounter_ += itr_;
156 totalItr_ += itr_;
157
158 return 0; // Success
159
160 } else{
161
162 // Not yet converged.
163 updateWFields();
164 timerMDE_.start();
165 system().compute();
166 ++mdeCounter_;
167 timerMDE_.stop();
168
169 }
170
171 }
172
173 // Failure: iteration counter itr reached maxItr without converging
174 timerTotal_.stop();
175
176 Log::file() << "Iterator failed to converge.\n";
177 return 1;
178
179 }
180
181 // Compute the residual for the current system state
182 template <int D>
184 {
185 const int nMonomer = system().mixture().nMonomer();
186
187 // Initialize resid to c field of species 0 minus 1
188 VecOp::subVS(resid_, system().c().rgrid(0), 1.0);
189
190 // Add other c fields to get SCF residual vector elements
191 for (int i = 1; i < nMonomer; i++) {
192 VecOp::addEqV(resid_, system().c().rgrid(i));
193 }
194 }
195
196 // update system w field using linear response approximation
197 template <int D>
198 void LrCompressor<D>::updateWFields(){
199 const int nMonomer = system().mixture().nMonomer();
200 const double vMonomer = system().mixture().vMonomer();
201
202 // Convert residual to Fourier Space
203 system().fft().forwardTransform(resid_, residK_);
204
205 // Compute change in fields using estimated Jacobian
206 VecOp::divEqVc(residK_, intraCorrelationK_, vMonomer);
207
208 // Convert back to real Space (destroys residK_)
209 system().fft().inverseTransformUnsafe(residK_, resid_);
210
211 // Update new fields
212 for (int i = 0; i < nMonomer; i++) {
213 VecOp::addVV(wFieldTmp_[i], system().w().rgrid(i), resid_);
214 }
215
216 // Set system r grid
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 return Reduce::maxAbs(a);
254 }
255
256 // Compute and return inner product of two vectors.
257 template <int D>
258 double LrCompressor<D>::dotProduct(RField<D> const & a,
259 RField<D> const & b)
260 {
261 UTIL_CHECK(a.capacity() == b.capacity());
262 return Reduce::innerProduct(a, b);
263 }
264
265 // Compute L2 norm of an RField
266 template <int D>
267 double LrCompressor<D>::norm(RField<D> const & a)
268 {
269 double normSq = dotProduct(a, a);
270 return sqrt(normSq);
271 }
272
273 // Compute and return the scalar error
274 template <int D>
275 double LrCompressor<D>::computeError(int verbose)
276 {
277 const int meshSize = system().domain().mesh().size();
278 double error = 0;
279
280 // Find max residual vector element
281 double maxRes = maxAbs(resid_);
282 // Find norm of residual vector
283 double normRes = norm(resid_);
284 // Find root-mean-squared residual element value
285 double rmsRes = normRes/sqrt(meshSize);
286 if (errorType_ == "maxResid") {
287 error = maxRes;
288 } else if (errorType_ == "normResid") {
289 error = normRes;
290 } else if (errorType_ == "rmsResid") {
291 error = rmsRes;
292 } else {
293 UTIL_THROW("Invalid iterator error type in parameter file.");
294 }
295
296 if (verbose > 1) {
297 Log::file() << "\n";
298 Log::file() << "Max Residual = " << Dbl(maxRes,15) << "\n";
299 Log::file() << "Residual Norm = " << Dbl(normRes,15) << "\n";
300 Log::file() << "RMS Residual = " << Dbl(rmsRes,15) << "\n";
301
302 // Check if calculation has diverged (normRes will be NaN)
303 UTIL_CHECK(!std::isnan(normRes));
304 }
305 return error;
306 }
307
308}
309}
310#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 setup()
Initialize just before entry to iterative loop.
void readParameters(std::istream &in)
Read all parameters and initialize.
void setClassName(const char *className)
Set class name string.
int compress()
Iterate pressure field to obtain partial saddle point.
void outputTimers(std::ostream &out)
Return compressor times contributions.
LrCompressor(System< D > &system)
Constructor.
Main class for calculations that represent one system.
Definition rpg/System.h:107
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 innerProduct(DeviceArray< cudaReal > const &a, DeviceArray< cudaReal > const &b)
Compute inner product of two real arrays (GPU kernel wrapper).
Definition Reduce.cu:891
cudaReal maxAbs(DeviceArray< cudaReal > const &in)
Get maximum absolute magnitude of array elements (GPU kernel wrapper).
Definition Reduce.cu:648
void addEqV(DeviceArray< cudaReal > &a, DeviceArray< cudaReal > const &b, const int beginIdA, const int beginIdB, const int n)
Vector addition in-place, a[i] += b[i], kernel wrapper (cudaReal).
Definition VecOp.cu:1632
void subVS(DeviceArray< cudaReal > &a, DeviceArray< cudaReal > const &b, cudaReal const c, const int beginIdA, const int beginIdB, const int n)
Vector subtraction, a[i] = b[i] - c, kernel wrapper (cudaReal).
Definition VecOp.cu:1304
void addVV(DeviceArray< cudaReal > &a, DeviceArray< cudaReal > const &b, DeviceArray< cudaReal > const &c, const int beginIdA, const int beginIdB, const int beginIdC, const int n)
Vector addition, a[i] = b[i] + c[i], kernel wrapper (cudaReal).
Definition VecOp.cu:1084
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
PSCF package top-level namespace.
Definition param_pc.dox:1
Utility classes for scientific computation.