PSCF v1.3.2
rpg/fts/compressor/AmCompressor.tpp
1#ifndef RPG_AM_COMPRESSOR_TPP
2#define RPG_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 "AmCompressor.h"
12#include <rpg/system/System.h>
13#include <rpg/solvers/Mixture.h>
14#include <rpg/field/Domain.h>
15#include <prdc/cuda/RField.h>
16#include <prdc/cuda/resources.h>
17#include <pscf/math/IntVec.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 isAllocated_(false)
32 { ParamComposite::setClassName("AmCompressor"); }
33
34 /*
35 * Destructor.
36 */
37 template <int D>
40
41 /*
42 * Read parameters from file.
43 */
44 template <int D>
45 void AmCompressor<D>::readParameters(std::istream& in)
46 {
47 // Default values
48 Base::maxItr_ = 100;
50 Base::errorType_ = "rms";
51 bool useLambdaRamp = false; // default
52
53 // Call parent class readParameters
56 Base::readMixingParameters(in, useLambdaRamp);
57 }
58
59 /*
60 * Initialize just before entry to iterative loop.
61 */
62 template <int D>
63 void AmCompressor<D>::setup(bool isContinuation)
64 {
65 const int nMonomer = system().mixture().nMonomer();
66 const int meshSize = system().domain().mesh().size();
67 const IntVec<D> dimensions = system().domain().mesh().dimensions();
68
69 // Allocate memory required by AM algorithm if not done earlier.
70 Base::setup(isContinuation);
71
72 // Allocate memory required by compressor if not done earlier.
73 if (!isAllocated_) {
74 w0_.allocate(nMonomer);
75 wFieldTmp_.allocate(nMonomer);
76 for (int i = 0; i < nMonomer; ++i) {
77 w0_[i].allocate(dimensions);
78 wFieldTmp_[i].allocate(dimensions);
79 }
80 isAllocated_ = true;
81 }
82
83 // Store value of initial guess chemical potential fields
84 for (int i = 0; i < nMonomer; ++i) {
85 VecOp::eqV(w0_[i], system().w().rgrid(i));
86 }
87 }
88
89 /*
90 * Main compressor function.
91 */
92 template <int D>
94 {
95 int solve = Base::solve();
96 return solve;
97 }
98
99 #if 0
100 /*
101 * Compute AM mixing parameter.
102 */
103 template<int D>
104 double AmCompressor<D>::computeLambda(double r)
105 { return 1.0; }
106 #endif
107
108 // Private virtual functions that interact with the parent System
109
110 /*
111 * Compute and return the number of elements in a field vector.
112 */
113 template <int D>
114 int AmCompressor<D>::nElements()
115 { return system().domain().mesh().size(); }
116
117 /*
118 * Does the system have an initial field guess?
119 */
120 template <int D>
121 bool AmCompressor<D>::hasInitialGuess()
122 { return system().w().hasData(); }
123
124 /*
125 * Get the current field from the system.
126 */
127 template <int D>
128 void AmCompressor<D>::getCurrent(VectorT& curr)
129 {
130 /*
131 * The field that we are adjusting is the Langrange multiplier field
132 * with number of grid pts components. The current value is the
133 * difference between w and w0_ for the first monomer (any monomer
134 * should give the same answer)
135 */
136 VecOp::subVV(curr, system().w().rgrid(0), w0_[0]);
137 }
138
139 /*
140 * Perform the main system computation (solve the MDE).
141 */
142 template <int D>
143 void AmCompressor<D>::evaluate()
144 {
145 system().compute();
147 }
148
149 /*
150 * Compute the residual for the current system state.
151 */
152 template <int D>
153 void AmCompressor<D>::getResidual(VectorT& resid)
154 {
155 const int nMonomer = system().mixture().nMonomer();
156
157 // Initialize resid to c field of species 0 minus 1
158 VecOp::subVS(resid, system().c().rgrid(0), 1.0);
159
160 // Add other c fields to get SCF residual vector elements
161 for (int i = 1; i < nMonomer; i++) {
162 VecOp::addEqV(resid, system().c().rgrid(i));
163 }
164 }
165
166 /*
167 * Update the current system field coordinates.
168 */
169 template <int D>
170 void AmCompressor<D>::update(VectorT& newGuess)
171 {
172 // Convert back to field format
173 const int nMonomer = system().mixture().nMonomer();
174
175 // New field is the w0_ + the newGuess for the Lagrange multiplier field
176 for (int i = 0; i < nMonomer; i++) {
177 VecOp::addVV(wFieldTmp_[i], w0_[i], newGuess);
178 }
179
180 // set system r grid
181 system().w().setRGrid(wFieldTmp_);
182 }
183
184 /*
185 * Do nothing output function.
186 */
187 template<int D>
188 void AmCompressor<D>::outputToLog()
189 {}
190
191 /*
192 * Write timing results to output stream
193 */
194 template<int D>
195 void AmCompressor<D>::outputTimers(std::ostream& out) const
196 {
197 out << "\n";
198 out << "Compressor times contributions:\n";
200 }
201
202 template<int D>
208
209 // Private virtual functions for vector math
210
211 /*
212 * Assign one array to another.
213 */
214 template <int D>
215 void AmCompressor<D>::setEqual(VectorT& a,
216 VectorT const & b)
217 {
218 UTIL_CHECK(b.capacity() == a.capacity());
219 VecOp::eqV(a, b);
220 }
221
222 /*
223 * Compute and return inner product of two vectors.
224 */
225 template <int D>
226 double AmCompressor<D>::dotProduct(VectorT const & a,
227 VectorT const & b)
228 {
229 UTIL_CHECK(a.capacity() == b.capacity());
230 return Reduce::innerProduct(a, b);
231 }
232
233 /*
234 * Compute and return maximum element of a vector.
235 */
236 template <int D>
237 double AmCompressor<D>::maxAbs(VectorT const & a)
238 { return Reduce::maxAbs(a); }
239
240 /*
241 * Compute the vector difference a = b - c
242 */
243 template <int D>
244 void AmCompressor<D>::subVV(VectorT& a,
245 VectorT const & b,
246 VectorT const & c)
247 { VecOp::subVV(a, b, c); }
248
249 /*
250 * Composite a += b*c for vectors a and b, scalar c
251 */
252 template <int D>
253 void AmCompressor<D>::addEqVc(VectorT& a,
254 VectorT const & b,
255 double c)
256 {
257 UTIL_CHECK(a.capacity() == b.capacity());
258 VecOp::addEqVc(a, b, c);
259 }
260
261}
262}
263#endif
virtual void setup(bool isContinuation)
void readMixingParameters(std::istream &in, bool useLambdaRamp=true)
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
int compress()
Compress to obtain partial saddle point w+.
AmCompressor(System< D > &system)
Constructor.
void setup(bool isContinuation)
Initialize just before entry to iterative loop.
void clearTimers()
Clear all timers (reset accumulated time to zero).
void readParameters(std::istream &in)
Read all parameters and initialize.
void outputTimers(std::ostream &out) const
Return compressor times contributions.
Base class for iterators that impose incompressibility.
int mdeCounter_
Count how many times MDE has been solved.
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 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.