PSCF v1.3
rpg/scft/sweep/Sweep.tpp
1#ifndef RPG_SWEEP_TPP
2/*
3* PSCF - Polymer Self-Consistent Field
4*
5* Copyright 2015 - 2025, The Regents of the University of Minnesota
6* Distributed under the terms of the GNU General Public License.
7*/
8
9#include "Sweep.h"
10#include <rpg/system/System.h>
11#include <rpg/scft/iterator/Iterator.h>
12#include <rpg/scft/ScftThermo.h>
13#include <prdc/environment/Environment.h>
14#include <pscf/inter/Interaction.h>
15#include <pscf/sweep/SweepTmpl.tpp>
16#include <util/misc/FileMaster.h>
17#include <util/misc/ioUtil.h>
18
19namespace Pscf {
20namespace Rpg {
21
22 using namespace Util;
23
24 // Maximum number of previous states = order of continuation + 1
25 #define RPG_HISTORY_CAPACITY 3
26
27 /*
28 * Default constructor.
29 */
30 template <int D>
32 : SweepTmpl< BasisFieldState<D> >(RPG_HISTORY_CAPACITY),
33 writeCRGrid_(false),
34 writeCBasis_(false),
35 writeWRGrid_(false),
36 systemPtr_(0)
37 {}
38
39 /*
40 * Constructor, creates association with parent system.
41 */
42 template <int D>
44 : SweepTmpl< BasisFieldState<D> >(RPG_HISTORY_CAPACITY),
45 writeCRGrid_(false),
46 writeCBasis_(false),
47 writeWRGrid_(false),
48 systemPtr_(&sys)
49 {
50 // Get specialized sweep parameters from Iterator
51 if (system().hasEnvironment()) {
52 addParameterTypes(system().environment().getParameterTypes());
53 }
54 }
55
56 /*
57 * Destructor.
58 */
59 template <int D>
62
63 /*
64 * Set association with a parent system.
65 */
66 template <int D>
68 { systemPtr_ = &system; }
69
70 /*
71 * Read parameters
72 */
73 template <int D>
74 void Sweep<D>::readParameters(std::istream& in)
75 {
76 // Call the base class's readParameters function.
78
79 // Read optional flags indicating which field types to output
80 readOptional(in, "writeCRGrid", writeCRGrid_);
81 readOptional(in, "writeCBasis", writeCBasis_);
82 readOptional(in, "writeWRGrid", writeWRGrid_);
83 }
84
85 /*
86 * Check allocation of one state object, allocate if necessary.
87 */
88 template <int D>
90 {
92 state.setSystem(system());
93 state.allocate();
94 state.unitCell() = system().domain().unitCell();
95 }
96
97 /*
98 * Setup operations at the beginning of a sweep.
99 */
100 template <int D>
102 {
103 initialize();
104 checkAllocation(trial_);
105
106 // Open log summary file
107 std::string fileName = baseFileName_;
108 fileName += "sweep.log";
109 system().fileMaster().openOutputFile(fileName, logFile_);
110 };
111
112 /*
113 * Set non-adjustable system parameters to new values.
114 *
115 * \param s path length coordinate, in range [0,1]
116 */
117 template <int D>
119 {
120 // Empty default implementation allows Sweep<D> to be compiled.
121 UTIL_THROW("Calling unimplemented function Sweep::setParameters");
122 };
123
124 /*
125 * Create guess for adjustable variables by polynomial extrapolation.
126 */
127 template <int D>
128 void Sweep<D>::extrapolate(double sNew)
129 {
130 UTIL_CHECK(historySize() > 0);
131
132 // If historySize() == 1, do nothing: Use previous system state
133 // as trial for the new state.
134
135 if (historySize() > 1) {
137
138 // Does the iterator allow a flexible unit cell ?
139 bool isFlexible = system().iterator().isFlexible();
140
141 // Compute coefficients of polynomial extrapolation to sNew
142 setCoefficients(sNew);
143
144 // Set extrapolated trial w fields
145 double coeff;
146 int nMonomer = system().mixture().nMonomer();
147 int nBasis = system().domain().basis().nBasis();
148 DArray<double>* newFieldPtr;
149 DArray<double>* oldFieldPtr;
150 int i, j, k;
151 for (i=0; i < nMonomer; ++i) {
152 newFieldPtr = &(trial_.field(i));
153
154 // Previous state k = 0 (most recent)
155 oldFieldPtr = &state(0).field(i);
156 coeff = c(0);
157 for (j=0; j < nBasis; ++j) {
158 (*newFieldPtr)[j] = coeff*(*oldFieldPtr)[j];
159 }
160
161 // Previous states k >= 1 (older)
162 for (k = 1; k < historySize(); ++k) {
163 oldFieldPtr = &state(k).field(i);
164 coeff = c(k);
165 for (j=0; j < nBasis; ++j) {
166 (*newFieldPtr)[j] += coeff*(*oldFieldPtr)[j];
167 }
168 }
169 }
170
171 // Make sure unitCellParameters_ is up to date with system
172 // (if we are sweeping in a lattice parameter, then the system
173 // parameters will be up-to-date but unitCellParameters_ wont be)
174 FSArray<double, 6> oldParameters = unitCellParameters_;
175 unitCellParameters_ = system().domain().unitCell().parameters();
176
177 // If isFlexible, then extrapolate the flexible unit cell parameters
178 if (isFlexible) {
179
180 double coeff;
181 double parameter;
182
183 #if 0
184 const FSArray<int,6> indices = system().iterator().flexibleParams();
185 const int nParameter = indices.size();
186
187 // Add contributions from all previous states
188 for (k = 0; k < historySize(); ++k) {
189 coeff = c(k);
190 for (int i = 0; i < nParameter; ++i) {
191 if (k == 0) {
192 unitCellParameters_[indices[i]] = 0;
193 }
194 parameter = state(k).unitCell().parameter(indices[i]);
195 unitCellParameters_[indices[i]] += coeff*parameter;
196 }
197 }
198 #endif
199
200 // Add contributions from all previous states
201 const int nParameter = unitCellParameters_.size();
202 for (k = 0; k < historySize(); ++k) {
203 coeff = c(k);
204 for (int i = 0; i < nParameter; ++i) {
205 if (k == 0) {
206 unitCellParameters_[i] = 0.0;
207 }
208 parameter = state(k).unitCell().parameter(i);
209 unitCellParameters_[i] += coeff*parameter;
210 }
211 }
212
213 }
214
215 // Reset trial_.unitCell() object
216 trial_.unitCell().setParameters(unitCellParameters_);
217
218 // Check if any unit cell parameters have changed
219 bool newCellParams(true);
220 for (int i = 0; i < oldParameters.size(); i++) {
221 if (fabs(oldParameters[i] - unitCellParameters_[i]) < 1e-10) {
222 newCellParams = false;
223 break;
224 }
225 }
226
227 // Transfer data from trial_ state to parent system
228 trial_.setSystemState(newCellParams);
229 }
230
231 };
232
233 /*
234 * Call current iterator to solve SCFT problem.
235 *
236 * Return 0 for sucessful solution, 1 on failure to converge.
237 */
238 template <int D>
239 int Sweep<D>::solve(bool isContinuation)
240 { return system().iterate(isContinuation); };
241
242 /*
243 * Reset system to previous solution after iterature failure.
244 *
245 * The implementation of this function should reset the system state
246 * to correspond to that stored in state(0).
247 */
248 template <int D>
250 {
251 bool isFlexible = system().iterator().isFlexible();
252 state(0).setSystemState(isFlexible);
253 }
254
255 /*
256 * Update state(0) and output data after successful convergence
257 *
258 * The implementation of this function should copy the current
259 * system state into state(0) and output any desired information
260 * about the current converged solution.
261 */
262 template <int D>
264 {
265 state(0).setSystem(system());
266 state(0).getSystemState();
267
268 // Output converged solution to several files
269 outputSolution();
270
271 // Output summary to log file
272 outputSummary(logFile_);
273
274 };
275
276 template <int D>
277 void Sweep<D>::outputSolution()
278 {
279 std::ofstream out;
280 std::string outFileName;
281 std::string indexString = toString(nAccept() - 1);
282
283 // Open parameter file, with thermodynamic properties at end
284 outFileName = baseFileName_;
285 outFileName += indexString;
286 outFileName += ".stt";
287 system().fileMaster().openOutputFile(outFileName, out);
288
289 // Write data file, with thermodynamic properties at end
290 system().writeParamNoSweep(out);
291 out << std::endl;
292 system().scft().write(out);
293 out.close();
294
295 // Write w fields
296 outFileName = baseFileName_;
297 outFileName += indexString;
298 outFileName += "_w";
299 outFileName += ".bf";
300 system().w().writeBasis(outFileName);
301
302 // Optionally write c rgrid files
303 if (writeCRGrid_) {
304 outFileName = baseFileName_;
305 outFileName += indexString;
306 outFileName += "_c";
307 outFileName += ".rf";
308 system().c().writeRGrid(outFileName);
309 }
310
311 // Optionally write c basis files
312 if (writeCBasis_) {
313 outFileName = baseFileName_;
314 outFileName += indexString;
315 outFileName += "_c";
316 outFileName += ".bf";
317 system().c().writeBasis(outFileName);
318 }
319
320 // Optionally write w rgrid files
321 if (writeWRGrid_) {
322 outFileName = baseFileName_;
323 outFileName += indexString;
324 outFileName += "_w";
325 outFileName += ".rf";
326 system().w().writeRGrid(outFileName);
327 }
328 }
329
330 template <int D>
331 void Sweep<D>::outputSummary(std::ostream& out)
332 {
333 int i = nAccept() - 1;
334 double sNew = s(0);
335 out << Int(i,5) << Dbl(sNew)
336 << Dbl(system().scft().fHelmholtz(),16)
337 << Dbl(system().scft().pressure(),16);
338 out << std::endl;
339 }
340
341 template <int D>
343 { logFile_.close(); }
344
345} // namespace Rpg
346} // namespace Pscf
347#endif
FieldState for fields in symmetry-adapted basis format.
virtual void getSolution()
Update state(0) and output data after successful convergence.
bool writeCBasis_
Whether to write concentration field files in basis format.
virtual int solve(bool isContinuation)
Call current iterator to solve SCFT problem.
virtual void readParameters(std::istream &in)
Read parameters from param file.
virtual void checkAllocation(BasisFieldState< D > &state)
Check allocation state of fields in one state, allocate if necessary.
virtual void extrapolate(double sNew)
Create a guess for adjustable variables by continuation.
ScalarParam< Type > & readOptional(std::istream &in, const char *label, Type &value)
Add and read a new optional ScalarParam < Type > object.
virtual void setup()
Setup operation at the beginning of a sweep.
bool writeWRGrid_
Whether to write real space potential field files.
bool writeCRGrid_
Whether to write real space concentration field files.
bool hasSystem()
Has an association with the parent System been set?
virtual void reset()
Reset system to previous solution after iterature failure.
virtual void setParameters(double sNew)=0
Set non-adjustable system parameters to new values.
Sweep()
Default Constructor.
System< D > & system()
Return the parent system by reference.
void setSystem(System< D > &system)
Set association with parent System.
virtual void cleanup()
Cleanup operation at the beginning of a sweep.
Main class, representing one complete system.
void addParameterTypes(GArray< ParameterType > paramTypes)
BasisFieldState< D > & state(int i)
Definition SweepTmpl.h:127
Dynamically allocatable contiguous array template.
Definition DArray.h:32
Wrapper for a double precision number, for formatted ostream output.
Definition Dbl.h:40
A fixed capacity (static) contiguous array with a variable logical size.
Definition FSArray.h:38
int size() const
Return logical size of this array (i.e., number of elements).
Definition FSArray.h:207
Wrapper for an int, for formatted ostream output.
Definition Int.h:37
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
Definition global.h:68
#define UTIL_THROW(msg)
Macro for throwing an Exception, reporting function, file and line number.
Definition global.h:49
std::string toString(int n)
Return string representation of an integer.
Definition ioUtil.cpp:52
SCFT and PS-FTS with real periodic fields (GPU)
PSCF package top-level namespace.
Definition param_pc.dox:1