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