PSCF v1.3.3
rpg/fts/compressor/LrAmCompressor.tpp
1#ifndef RPG_LR_AM_COMPRESSOR_TPP
2#define RPG_LR_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 <rpg/solvers/Mixture.h>
14#include <rpg/field/Domain.h>
15#include <prdc/cuda/FFT.h>
16#include <prdc/cuda/resources.h>
17#include <pscf/mesh/MeshIterator.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 intra_(system),
32 isIntraCalculated_(false),
33 isAllocated_(false)
34 { ParamComposite::setClassName("LrAmCompressor"); }
35
36 /*
37 * Destructor.
38 */
39 template <int D>
42
43 /*
44 * Read parameters from file.
45 */
46 template <int D>
47 void LrAmCompressor<D>::readParameters(std::istream& in)
48 {
49 // Default values
50 Base::maxItr_ = 100;
52 Base::errorType_ = "rms";
53
54 // Call base class read methods
57 }
58
59 /*
60 * Initialize just before entry to iterative loop.
61 */
62 template <int D>
63 void LrAmCompressor<D>::setup(bool isContinuation)
64 {
65 const int nMonomer = system().mixture().nMonomer();
66 const int meshSize = system().domain().mesh().size();
67 IntVec<D> const & dimensions = system().domain().mesh().dimensions();
68
69 // Allocate memory required by AM algorithm if not done earlier.
70 Base::setup(isContinuation);
71
72 FFT<D>::computeKMesh(dimensions, kMeshDimensions_, kSize_);
73
74 // Compute Fourier space kMeshDimensions_
75 for (int i = 0; i < D; ++i) {
76 if (i < D - 1) {
77 kMeshDimensions_[i] = dimensions[i];
78 } else {
79 kMeshDimensions_[i] = dimensions[i]/2 + 1;
80 }
81 }
82
83 // Compute number of points in k-space grid
84 kSize_ = 1;
85 for (int i = 0; i < D; ++i) {
86 kSize_ *= kMeshDimensions_[i];
87 }
88
89 // Allocate memory required by compressor if not done earlier.
90 if (!isAllocated_) {
91 w0_.allocate(nMonomer);
92 wFieldTmp_.allocate(nMonomer);
93 resid_.allocate(dimensions);
94 residK_.allocate(dimensions);
95 intraCorrelationK_.allocate(kMeshDimensions_);
96 for (int i = 0; i < nMonomer; ++i) {
97 w0_[i].allocate(dimensions);
98 wFieldTmp_[i].allocate(dimensions);
99 }
100
101 isAllocated_ = true;
102 }
103
104 // Compute intraCorrelation
105 if (!isIntraCalculated_){
106 intra_.computeIntraCorrelations(intraCorrelationK_);
107 isIntraCalculated_ = true;
108 }
109
110 // Store value of initial guess chemical potential fields
111 for (int i = 0; i < nMonomer; ++i) {
112 VecOp::eqV(w0_[i], system().w().rgrid(i));
113 }
114
115 }
116
117 /*
118 * Main function - iteratively adjust the pressure field.
119 */
120 template <int D>
122 {
123 int solve = Base::solve();
124 return solve;
125 }
126
127 // Protected virtual function for AM algorithm operation
128
129 template <int D>
130 void
131 LrAmCompressor<D>::addCorrection(
132 VectorT& fieldTrial,
133 VectorT const & resTrial)
134 {
135 int n = fieldTrial.capacity();
136 const double vMonomer = system().mixture().vMonomer();
137
138 // Convert resTrial to RField<D> type
139 VecOp::eqV(resid_, resTrial);
140
141 // Convert residual to Fourier Space
142 system().domain().fft().forwardTransform(resid_, residK_);
143
144 // Combine with Linear response factor to update second step
145 VecOp::divEqVc(residK_, intraCorrelationK_, vMonomer);
146
147 // Convert back to real space (destroys residK_)
148 system().domain().fft().inverseTransformUnsafe(residK_, resid_);
149
150 // fieldTrial += resid_
151 VecOp::addEqVc(fieldTrial, resid_, 1.0);
152 }
153
154 // Private virtual functions that interact with parent System
155
156 /*
157 * Get the number of elements in a field vector.
158 */
159 template <int D>
160 int LrAmCompressor<D>::nElements()
161 { return system().domain().mesh().size(); }
162
163 /*
164 * Does the system have an initial field guess?
165 */
166 template <int D>
167 bool LrAmCompressor<D>::hasInitialGuess()
168 { return system().w().hasData(); }
169
170 /*
171 * Get the current field from the system.
172 */
173 template <int D>
174 void LrAmCompressor<D>::getCurrent(VectorT& curr)
175 {
176 /*
177 * The field that we are adjusting is the Langrange multiplier field
178 * with number of grid pts components.The current value is the difference
179 * between w and w0_ for the first monomer (any monomer should give the
180 * same answer)
181 */
182 VecOp::subVV(curr, system().w().rgrid(0), w0_[0]);
183 }
184
185 /*
186 * Perform the main system computation (solve the MDE).
187 */
188 template <int D>
189 void LrAmCompressor<D>::evaluate()
190 {
191 system().compute();
193 }
194
195 /*
196 * Compute the residual vector for the current system state.
197 */
198 template <int D>
199 void LrAmCompressor<D>::getResidual(VectorT& resid)
200 {
201 const int nMonomer = system().mixture().nMonomer();
202
203 // Initialize resid to c field of species 0 minus 1
204 VecOp::subVS(resid, system().c().rgrid(0), 1.0);
205
206 // Add other c fields to get SCF residual vector elements
207 for (int i = 1; i < nMonomer; i++) {
208 VecOp::addEqV(resid, system().c().rgrid(i));
209 }
210 }
211
212 /*
213 * Update the system state to reflect new state vector.
214 */
215 template <int D>
216 void LrAmCompressor<D>::update(VectorT& newGuess)
217 {
218 // Convert back to field format
219 const int nMonomer = system().mixture().nMonomer();
220
221 // New field is w0_[i] + newGuess for the Lagrange multiplier field
222 for (int i = 0; i < nMonomer; i++) {
223 VecOp::addVV(wFieldTmp_[i], w0_[i], newGuess);
224 }
225
226 // Set system r grid
227 system().w().setRGrid(wFieldTmp_);
228 }
229
230 /*
231 * Do-nothing output function.
232 */
233 template<int D>
234 void LrAmCompressor<D>::outputToLog()
235 {}
236
237 /*
238 * Output timer results.
239 */
240 template<int D>
241 void LrAmCompressor<D>::outputTimers(std::ostream& out) const
242 {
243 out << "\n";
244 out << "LrAmCompressor time contributions:\n";
246 }
247
248 /*
249 * Clear timers and MDE counter.
250 */
251 template<int D>
257
258 // Private virtual functions for vector math
259
260 /*
261 * Vector assignment, a = b.
262 */
263 template <int D>
264 void LrAmCompressor<D>::setEqual(VectorT& a,
265 VectorT const & b)
266 {
267 UTIL_CHECK(b.capacity() == a.capacity());
268 VecOp::eqV(a, b);
269 }
270
271 /*
272 * Compute and return inner product of two vectors.
273 */
274 template <int D>
275 double LrAmCompressor<D>::dotProduct(VectorT const & a,
276 VectorT const & b)
277 {
278 UTIL_CHECK(a.capacity() == b.capacity());
279 return Reduce::innerProduct(a, b);
280 }
281
282 /*
283 * Compute and return maximum element of a vector.
284 */
285 template <int D>
286 double LrAmCompressor<D>::maxAbs(VectorT const & a)
287 { return Reduce::maxAbs(a); }
288
289 /*
290 * Compute the vector difference a = b - c
291 */
292 template <int D>
293 void LrAmCompressor<D>::subVV(VectorT& a,
294 VectorT const & b,
295 VectorT const & c)
296 { VecOp::subVV(a, b, c); }
297
298 /*
299 * Compute a += b*c for vectors a and b, scalar c
300 */
301 template <int D>
302 void LrAmCompressor<D>::addEqVc(VectorT& a,
303 VectorT const & b,
304 double c)
305 { VecOp::addEqVc(a, b, c); }
306
307}
308}
309#endif
virtual void setup(bool isContinuation)
void outputTimers(std::ostream &out) const
virtual void readParameters(std::istream &in)
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.
int mdeCounter_
Count how many times MDE has been solved.
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 a complete physical system.
void setClassName(const char *className)
Set class name string.
File containing preprocessor macros for error handling.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
Definition global.h:68
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 subVV(Array< double > &a, Array< double > const &b, Array< double > const &c)
Vector subtraction, a[i] = b[i] - c[i].
Definition VecOp.cpp:81
void addEqV(Array< double > &a, Array< double > const &b)
Vector addition in-place, a[i] += b[i].
Definition VecOp.cpp:199
void eqV(Array< double > &a, Array< double > const &b)
Vector assignment, a[i] = b[i].
Definition VecOp.cpp:23
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
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.