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