PSCF v1.2
rpg/fts/compressor/LrAmCompressor.tpp
1#ifndef RPG_LR_POST_AM_COMPRESSOR_TPP
2#define RPG_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 <rpg/System.h>
13#include <rpg/fts/compressor/intra/IntraCorrelation.h>
14#include <prdc/cuda/resources.h>
15#include <pscf/mesh/MeshIterator.h>
16#include <util/global.h>
17
18namespace Pscf {
19namespace Rpg{
20
21 using namespace Util;
22
23 // Constructor
24 template <int D>
26 : Compressor<D>(system),
27 intra_(system),
28 isIntraCalculated_(false),
29 isAllocated_(false)
30 { setClassName("LrAmCompressor"); }
31
32 // Destructor
33 template <int D>
36
37 // Read parameters from file
38 template <int D>
39 void LrAmCompressor<D>::readParameters(std::istream& in)
40 {
41 // Call parent class readParameters
44 }
45
46 // Initialize just before entry to iterative loop.
47 template <int D>
48 void LrAmCompressor<D>::setup(bool isContinuation)
49 {
50 const int nMonomer = system().mixture().nMonomer();
51 const int meshSize = system().domain().mesh().size();
52 IntVec<D> const & dimensions = system().mesh().dimensions();
53
54 // Allocate memory required by AM algorithm if not done earlier.
56
57 // Compute Fourier space kMeshDimensions_
58 for (int i = 0; i < D; ++i) {
59 if (i < D - 1) {
60 kMeshDimensions_[i] = dimensions[i];
61 } else {
62 kMeshDimensions_[i] = dimensions[i]/2 + 1;
63 }
64 }
65
66 // Compute number of points in k-space grid
67 kSize_ = 1;
68 for (int i = 0; i < D; ++i) {
69 kSize_ *= kMeshDimensions_[i];
70 }
71
72 // Allocate memory required by compressor if not done earlier.
73 if (!isAllocated_) {
74 w0_.allocate(nMonomer);
75 wFieldTmp_.allocate(nMonomer);
76 resid_.allocate(dimensions);
77 residK_.allocate(dimensions);
78 intraCorrelationK_.allocate(kMeshDimensions_);
79 for (int i = 0; i < nMonomer; ++i) {
80 w0_[i].allocate(dimensions);
81 wFieldTmp_[i].allocate(dimensions);
82 }
83
84 isAllocated_ = true;
85 }
86
87 // Compute intraCorrelation
88 if (!isIntraCalculated_){
89 intra_.computeIntraCorrelations(intraCorrelationK_);
90 isIntraCalculated_ = true;
91 }
92
93 // Store value of initial guess chemical potential fields
94 for (int i = 0; i < nMonomer; ++i) {
95 VecOp::eqV(w0_[i], system().w().rgrid(i));
96 }
97
98 }
99
100 template <int D>
102 {
104 //mdeCounter_ = AmIteratorTmpl<Compressor<D>,DArray<double>>::totalItr();
105 return solve;
106 }
107
108 // Assign one array to another
109 template <int D>
111 DeviceArray<cudaReal> const & b)
112 {
113 UTIL_CHECK(b.capacity() == a.capacity());
114 VecOp::eqV(a, b);
115 }
116
117 // Compute and return inner product of two vectors.
118 template <int D>
119 double LrAmCompressor<D>::dotProduct(DeviceArray<cudaReal> const & a,
120 DeviceArray<cudaReal> const & b)
121 {
122 UTIL_CHECK(a.capacity() == b.capacity());
123 return Reduce::innerProduct(a, b);
124 }
125
126 // Compute and return maximum element of a vector.
127 template <int D>
128 double LrAmCompressor<D>::maxAbs(DeviceArray<cudaReal> const & a)
129 {
130 return Reduce::maxAbs(a);
131 }
132
133 // Update basis
134 template <int D>
135 void
136 LrAmCompressor<D>::updateBasis(RingBuffer< DeviceArray<cudaReal> > & basis,
137 RingBuffer< DeviceArray<cudaReal> > const & hists)
138 {
139 // Make sure at least two histories are stored
140 UTIL_CHECK(hists.size() >= 2);
141
142 basis.advance();
143 if (basis[0].isAllocated()) {
144 UTIL_CHECK(basis[0].capacity() == hists[0].capacity());
145 } else {
146 basis[0].allocate(hists[0].capacity());
147 }
148
149 VecOp::subVV(basis[0], hists[0], hists[1]);
150 }
151
152 template <int D>
153 void
154 LrAmCompressor<D>::addHistories(DeviceArray<cudaReal>& trial,
155 RingBuffer<DeviceArray<cudaReal> > const & basis,
156 DArray<double> coeffs,
157 int nHist)
158 {
159 for (int i = 0; i < nHist; i++) {
160 VecOp::addEqVc(trial, basis[i], -1.0 * coeffs[i]);
161 }
162 }
163
164 template<int D>
165 double LrAmCompressor<D>::computeLambda(double r)
166 {
167 return 1.0;
168 }
169
170 template <int D>
171 void LrAmCompressor<D>::addPredictedError(DeviceArray<cudaReal>& fieldTrial,
172 DeviceArray<cudaReal> const & resTrial,
173 double lambda)
174 {
175 int n = fieldTrial.capacity();
176 const double vMonomer = system().mixture().vMonomer();
177
178 // Convert resTrial to RField<D> type
179 VecOp::eqV(resid_, resTrial);
180
181 // Convert residual to Fourier Space
182 system().fft().forwardTransform(resid_, residK_);
183
184 // Combine with Linear response factor to update second step
185 VecOp::divEqVc(residK_, intraCorrelationK_, vMonomer);
186
187 // Convert back to real space (destroys residK_)
188 system().fft().inverseTransformUnsafe(residK_, resid_);
189
190 // fieldTrial += resid_ * lambda
191 VecOp::addEqVc(fieldTrial, resid_, lambda);
192 }
193
194 // Does the system have an initial field guess?
195 template <int D>
196 bool LrAmCompressor<D>::hasInitialGuess()
197 { return system().w().hasData(); }
198
199 // Compute and return the number of elements in a field vector
200 template <int D>
201 int LrAmCompressor<D>::nElements()
202 { return system().domain().mesh().size(); }
203
204 // Get the current field from the system
205 template <int D>
206 void LrAmCompressor<D>::getCurrent(DeviceArray<cudaReal>& curr)
207 {
208 /*
209 * The field that we are adjusting is the Langrange multiplier field
210 * with number of grid pts components.The current value is the difference
211 * between w and w0_ for the first monomer (any monomer should give the
212 * same answer)
213 */
214 VecOp::subVV(curr, system().w().rgrid(0), w0_[0]);
215 }
216
217 // Perform the main system computation (solve the MDE)
218 template <int D>
219 void LrAmCompressor<D>::evaluate()
220 {
221 system().compute();
222 ++mdeCounter_;
223 }
224
225 // Compute the residual for the current system state
226 template <int D>
227 void LrAmCompressor<D>::getResidual(DeviceArray<cudaReal>& resid)
228 {
229 const int nMonomer = system().mixture().nMonomer();
230
231 // Initialize resid to c field of species 0 minus 1
232 VecOp::subVS(resid, system().c().rgrid(0), 1.0);
233
234 // Add other c fields to get SCF residual vector elements
235 for (int i = 1; i < nMonomer; i++) {
236 VecOp::addEqV(resid, system().c().rgrid(i));
237 }
238 }
239
240 // Update the current system field coordinates
241 template <int D>
242 void LrAmCompressor<D>::update(DeviceArray<cudaReal>& newGuess)
243 {
244 // Convert back to field format
245 const int nMonomer = system().mixture().nMonomer();
246
247 // New field is the w0_ + the newGuess for the Lagrange multiplier field
248 for (int i = 0; i < nMonomer; i++) {
249 VecOp::addVV(wFieldTmp_[i], w0_[i], newGuess);
250 }
251
252 // Set system r grid
253 system().setWRGrid(wFieldTmp_);
254 }
255
256 template<int D>
257 void LrAmCompressor<D>::outputToLog()
258 {}
259
260 template<int D>
261 void LrAmCompressor<D>::outputTimers(std::ostream& out)
262 {
263 // Output timing results, if requested.
264 out << "\n";
265 out << "LrAmCompressor time contributions:\n";
267 }
268
269 // Clear timers and MDE counter
270 template<int D>
272 {
274 mdeCounter_ = 0;
275 }
276
277}
278}
279#endif
Template for Anderson mixing iterator algorithm.
Dynamic array on the GPU device with aligned data.
Definition rpg/System.h:32
int capacity() const
Return allocated capacity.
An IntVec<D, T> is a D-component vector of elements of integer type T.
Definition IntVec.h:27
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.
int compress()
Compress to obtain partial saddle point w+.
void clearTimers()
Clear all timers (reset accumulated time to zero).
void setup(bool isContinuation)
Initialize just before entry to iterative loop.
LrAmCompressor(System< D > &system)
Constructor.
void outputTimers(std::ostream &out)
Return compressor times contributions.
Main class for calculations that represent one system.
Definition rpg/System.h:107
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 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 eqV(DeviceArray< cudaReal > &a, DeviceArray< cudaReal > const &b, const int beginIdA, const int beginIdB, const int n)
Vector assignment, a[i] = b[i], kernel wrapper (cudaReal).
Definition VecOp.cu:1020
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 subVV(DeviceArray< cudaReal > &a, DeviceArray< cudaReal > const &b, DeviceArray< cudaReal > const &c, const int beginIdA, const int beginIdB, const int beginIdC, const int n)
Vector subtraction, a[i] = b[i] - c[i], kernel wrapper (cudaReal).
Definition VecOp.cu:1228
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
void addEqVc(DeviceArray< cudaReal > &a, DeviceArray< cudaReal > const &b, cudaReal const c)
Vector addition in-place w/ coefficient, a[i] += b[i] * c, kernel wrapper.
Definition VecOpMisc.cu:343
PSCF package top-level namespace.
Definition param_pc.dox:1
Utility classes for scientific computation.