PSCF v1.3
PolymerSpecies.cpp
1/*
2* PSCF - Polymer Self-Consistent Field
3*
4* Copyright 2015 - 2025, The Regents of the University of Minnesota
5* Distributed under the terms of the GNU General Public License.
6*/
7
8#include "PolymerSpecies.h"
9
10#include <pscf/chem/Edge.h>
11#include <util/containers/GArray.h>
12#include <util/containers/FArray.h>
13#include <util/containers/DMatrix.h>
14
15#include <cmath>
16
17namespace Pscf
18{
19
20 /*
21 * Constructor.
22 */
24 : Species(),
25 vertices_(),
26 propagatorIds_(),
27 paths_(),
28 nBlock_(0),
29 nVertex_(0),
30 nPropagator_(0),
31 type_(PolymerType::Linear)
32 { setClassName("PolymerSpecies"); }
33
34 /*
35 * Destructor.
36 */
39
40 /*
41 * Read parameter file block.
42 */
43 void PolymerSpecies::readParameters(std::istream& in)
44 {
45 // Read polymer type (linear by default)
46 type_ = PolymerType::Linear;
47 readOptional<PolymerType::Enum>(in, "type", type_);
48
49 read<int>(in, "nBlock", nBlock_);
50
51 // Note: For any acyclic graph, nVertex = nBlock + 1
52 nVertex_ = nBlock_ + 1;
53
54 // Allocate all arrays
56 vertices_.allocate(nVertex_);
57 propagatorIds_.allocate(2*nBlock_);
58
59 // Set block id and polymerType for all blocks
60 for (int edgeId = 0; edgeId < nBlock_; ++edgeId) {
61 edge(edgeId).setId(edgeId);
62 edge(edgeId).setPolymerType(type_);
63 }
64
65 // Set sequential ids for all vertices
66 for (int vertexId = 0; vertexId < nVertex_; ++vertexId) {
67 vertices_[vertexId].setId(vertexId);
68 }
69
70 // If this is a linear polymer, set all block vertex Ids:
71 if (type_ == PolymerType::Linear) {
72 // In a linear chain, block i connects vertices i and i+1.
73 for (int blockId = 0; blockId < nBlock_; ++blockId) {
74 edge(blockId).setVertexIds(blockId, blockId + 1);
75 }
76 }
77
78 // Read array of block data from parameter file
79 readBlocks(in);
80
81 // Read phi or mu (but not both) and set ensemble accordingly
83
84 // Reading of parameter file is now complete
85
86 // Add edges to attached vertices
87 int vertexId0, vertexId1;
88 Edge* edgePtr;
89 for (int edgeId = 0; edgeId < nBlock_; ++edgeId) {
90 edgePtr = &(edge(edgeId));
91 vertexId0 = edgePtr->vertexId(0);
92 vertexId1 = edgePtr->vertexId(1);
93 vertices_[vertexId0].addEdge(*edgePtr);
94 vertices_[vertexId1].addEdge(*edgePtr);
95 }
96
97 // Polymer graph topology is now fully specified.
98
99 // Construct a plan for the order in which block propagators
100 // should be computed when solving the MDE.
101 makePlan();
102
103 // Construct vertex-to-vertex path signposts
104 makePaths();
105
106 // Check internal consistency of edge and vertex data
107 isValid();
108
109 }
110
111 /*
112 * Make a plan for the order of computation of block propagators.
113 */
115 {
116 if (nPropagator_ != 0) {
117 UTIL_THROW("nPropagator !=0 on entry");
118 }
119
120 // Allocate and initialize isFinished matrix
121 DMatrix<bool> isFinished;
122 isFinished.allocate(nBlock_, 2);
123 for (int iBlock = 0; iBlock < nBlock_; ++iBlock) {
124 for (int iDirection = 0; iDirection < 2; ++iDirection) {
125 isFinished(iBlock, iDirection) = false;
126 }
127 }
128
129 Pair<int> propId;
130 Vertex* inVertexPtr = 0;
131 int inVertexId = -1;
132 bool isReady;
133 while (nPropagator_ < nBlock_*2) {
134 for (int iBlock = 0; iBlock < nBlock_; ++iBlock) {
135 for (int iDirection = 0; iDirection < 2; ++iDirection) {
136 if (isFinished(iBlock, iDirection) == false) {
137 inVertexId = edge(iBlock).vertexId(iDirection);
138 inVertexPtr = &vertices_[inVertexId];
139 isReady = true;
140 for (int j = 0; j < inVertexPtr->size(); ++j) {
141 propId = inVertexPtr->inPropagatorId(j);
142 if (propId[0] != iBlock) {
143 if (!isFinished(propId[0], propId[1])){
144 isReady = false;
145 break;
146 }
147 }
148 }
149 if (isReady) {
150 propagatorIds_[nPropagator_][0] = iBlock;
151 propagatorIds_[nPropagator_][1] = iDirection;
152 isFinished(iBlock, iDirection) = true;
153 ++nPropagator_;
154 }
155 }
156 }
157 }
158 }
159
161
162 /*
163 * Identify paths between vertices.
164 */
166 {
167 UTIL_CHECK(nVertex_ > 0);
168
169 // Local variables
170 DArray< GArray< FArray<int, 3> > > sendPaths;
171 DArray< GArray< FArray<int, 3> > > recvPaths;
174 Pair<int> pId;
175 int is; // sender vertex id
176 int ir; // receiver vertex id
177 int ib; // bond id
178 int io; // outgoing (source) direction id
179 int ii; // incoming (receiver) direction id
180 int j, k, n;
181
182 // Allocate and clear path containers
183 sendPaths.allocate(nVertex_);
184 recvPaths.allocate(nVertex_);
185 oldPaths.allocate(nVertex_);
186 for (is = 0; is < nVertex_; ++is) {
187 sendPaths[is].clear();
188 recvPaths[is].clear();
189 oldPaths[is].clear();
190 }
191
192 // Initialize sendPaths with path-to-self entries
193 for (is = 0; is < nVertex_; ++is) {
194 path[0] = is;
195 path[1] = -1;
196 path[2] = -1;
197 sendPaths[is].append(path);
198 }
199
200 // While loop to completion
201 bool done = false;
202 int iter = 0;
203 while (!done) {
204 UTIL_CHECK(iter < nBlock_ + 2);
205
206 // Check that recvPaths container is empty
207 for (ir = 0; is < nVertex_; ++is) {
208 UTIL_CHECK(recvPaths[ir].size() == 0);
209 }
210
211 // Loop over source vertices for sending
212 done = true;
213 for (is = 0; is < nVertex_; ++is) {
214
215 // Send any sendPaths data to neighbors
216 int n = sendPaths[is].size();
217 if (n > 0) {
218 done = false;
219 Vertex const & sender = vertex(is);
220
221 // Send sendPaths data to all new neighbors
222 for (j = 0; j < sender.size(); ++j) {
223 pId = sender.outPropagatorId(j);
224 ib = pId[0]; // block identifier
225 io = pId[1]; // outgoing direction identifier
226 if (io == 0) {
227 UTIL_CHECK(edge(ib).vertexId(0) == is);
228 ii = 1;
229 } else {
230 UTIL_CHECK(edge(ib).vertexId(1) == is);
231 ii = 0;
232 }
233 ir = edge(ib).vertexId(ii);
234 for (k = 0; k < n; ++k) {
235 path = sendPaths[is][k];
236 // Send unless just received along same bond
237 if (ib != path[1]) {
238 path[1] = ib;
239 path[2] = ii;
240 recvPaths[ir].append(path);
241 }
242 }
243 }
244
245 } // if (n > 0)
246 } // Loop over source vertices
247
248 // Loop over vertices to transfer data structures
249 for (is = 0; is < nVertex_; ++is) {
250
251 // Move sendPaths to oldPaths, clear sendPaths
252 n = sendPaths[is].size();
253 if (n > 0) {
254 for (k = 0; k < n; ++k) {
255 path = sendPaths[is][k];
256 oldPaths[is].append(path);
257 }
258 }
259 sendPaths[is].clear();
260
261 // Move recvPaths to sendPaths, clear recvPaths
262 n = recvPaths[is].size();
263 if (n > 0) {
264 for (k = 0; k < n; ++k) {
265 sendPaths[is].append(recvPaths[is][k]);
266 }
267 }
268 recvPaths[is].clear();
269
270 }
271
272 ++iter;
273 } // while not done
274
275 // Allocate and initialize member variable paths_
276 paths_.allocate(nVertex_);
277 for (is = 0; is < nVertex_; ++is) {
278 paths_[is].allocate(nVertex_);
279 for (ir = 0; ir < nVertex_; ++ir) {
280 paths_[is][ir][0] = -1;
281 paths_[is][ir][1] = -1;
282 }
283 }
284
285 // Assign values to all elements of paths_ container
286 for (is = 0; is < nVertex_; ++is) {
287 n = oldPaths[is].size();
288 UTIL_CHECK(n == nVertex_);
289 for (k = 0; k < n; ++k) {
290 path = oldPaths[is][k];
291 ir = path[0]; // id of target vertex
292 UTIL_CHECK(ir >= 0);
293 UTIL_CHECK(ir < nVertex_);
294 // Check that element was not set previously
295 UTIL_CHECK(paths_[is][ir][0] == -1);
296 UTIL_CHECK(paths_[is][ir][1] == -1);
297 //std::cout << std::endl << is << " " << ir << " "
298 // << path[1] << " " << path[2];
299 if (ir == is) {
300 UTIL_CHECK(path[1] == -1);
301 UTIL_CHECK(path[2] == -1);
302 } else {
303 UTIL_CHECK(path[1] >= 0);
304 UTIL_CHECK(path[1] < nBlock_);
305 UTIL_CHECK(path[2] >= 0);
306 UTIL_CHECK(path[2] < 2);
307 }
308 paths_[is][ir][0] = path[1];
309 paths_[is][ir][1] = path[2];
310 }
311 }
312
313 #if 0
314 std::cout << std::endl << "Paths:" ;
315 for (is = 0; is < nVertex_; ++is) {
316 for (ir = 0; ir < nVertex_; ++ir) {
317 pId = paths_[is][ir];
318 std::cout << std::endl;
319 std::cout << is << ir << pId[0] << pId[1];
320 }
321 }
322 std::cout << std::endl;
323 #endif
324
325 }
326
327 /*
328 * Check consistency of internal data structures.
329 */
330 void PolymerSpecies::isValid()
331 {
332 Pair<int> pair;
333 int ib, iv, ip, id, iv0, iv1, n;
334
335 // Check validity of ids owned by edges
336 for (ib = 0; ib < nBlock_; ++ib) {
337 UTIL_CHECK(edge(ib).id() == ib);
338 iv0 = edge(ib).vertexId(0);
339 iv1 = edge(ib).vertexId(1);
340 UTIL_CHECK(iv0 != iv1);
341 UTIL_CHECK(iv0 >= 0);
342 UTIL_CHECK(iv0 < nVertex_);
343 UTIL_CHECK(iv1 >= 0);
344 UTIL_CHECK(iv1 < nVertex_);
345 }
346
347 // Check consistency of vertex::outPropagatorId
348 for (iv = 0; iv < nVertex_; ++iv) {
349 UTIL_CHECK(vertex(iv).id() == iv);
350 n = vertex(iv).size();
351 for (ip = 0; ip < n; ++ip) {
352 pair = vertex(iv).outPropagatorId(ip);
353 ib = pair[0];
354 id = pair[1];
355 UTIL_CHECK(ib >= 0);
356 UTIL_CHECK(ib < nBlock_);
357 UTIL_CHECK(id >= 0);
358 UTIL_CHECK(id < 2);
359 UTIL_CHECK(edge(ib).vertexId(id) == iv);
360 }
361 }
362
363 // Check consistency of vertex::inPropagatorId
364 for (iv = 0; iv < nVertex_; ++iv) {
365 UTIL_CHECK(vertex(iv).id() == iv);
366 n = vertex(iv).size();
367 for (ip = 0; ip < n; ++ip) {
368 pair = vertex(iv).inPropagatorId(ip);
369 ib = pair[0];
370 id = pair[1];
371 UTIL_CHECK(ib >= 0);
372 UTIL_CHECK(ib < nBlock_);
373 UTIL_CHECK(id >= 0);
374 UTIL_CHECK(id < 2);
375 if (id == 0) {
376 UTIL_CHECK(edge(ib).vertexId(1) == iv);
377 } else {
378 UTIL_CHECK(edge(ib).vertexId(0) == iv);
379 }
380 }
381 }
382
383 // Check consistency of vertex ids in paths
384 for (iv0 = 0; iv0 < nVertex_; ++iv0) {
385 for (iv1 = 0; iv1 < nVertex_; ++iv1) {
386 pair = path(iv0, iv1);
387 ib = pair[0];
388 id = pair[1];
389 if (iv0 == iv1) {
390 UTIL_CHECK(ib == -1);
391 UTIL_CHECK(id == -1);
392 } else {
393 UTIL_CHECK(edge(ib).vertexId(id) == iv0);
394 }
395 }
396 }
397
398 }
399
400 /*
401 * Total length of all blocks = volume / reference volume
402 */
404 {
406 double value = 0.0;
407 for (int blockId = 0; blockId < nBlock_; ++blockId) {
408 value += edge(blockId).length();
409 }
410 return value;
411 }
412
413 /*
414 * Total number of beads in all blocks.
415 */
417 {
419 int value = 0.0;
420 for (int blockId = 0; blockId < nBlock_; ++blockId) {
421 value += edge(blockId).nBead();
422 }
423 return value;
424 }
425
426}
Descriptor for a block within a block polymer.
Definition Edge.h:59
void setPolymerType(PolymerType::Enum type)
Set the type of the parent polymer (branched or linear).
Definition Edge.cpp:77
void setVertexIds(int vertexId0, int vertexId1)
Set indices of associated vertices.
Definition Edge.cpp:50
int nBead() const
Get the number of beads in this block, in the bead model.
Definition Edge.h:302
int vertexId(int i) const
Get the id of one associated vertex.
Definition Edge.h:296
void setId(int id)
Set the id for this block.
Definition Edge.cpp:38
double length() const
Get the length of this block, in the thread model.
Definition Edge.h:311
const Vertex & vertex(int id) const
Get a specified Vertex by const reference.
virtual void readBlocks(std::istream &in)=0
Read array of blocks from parameter file.
PolymerSpecies()
Constructor.
void makePaths()
Create a matrix of vertex-to-vertex path signposts.
Pair< int > const & path(int is, int it) const
Get an id for a propagator from one vertex towards a target.
virtual void allocateBlocks()=0
Allocate array of blocks.
virtual Edge & edge(int id)=0
Get a specified Edge (block descriptor) by non-const reference.
virtual void readParameters(std::istream &in)
Read parameters and initialize.
virtual void makePlan()
Make a plan for order in which propagators should be computed.
~PolymerSpecies()
Destructor.
int nBead() const
Total number of beads in the polymer (bead model).
double length() const
Sum of the lengths of all blocks in the polymer (thread model).
virtual void readParameters(std::istream &in)
Read phi or mu (but not both) and set ensemble accordingly.
Definition Species.cpp:28
Species()
Default constructor.
Definition Species.cpp:18
A junction or chain end in a block polymer.
Definition Vertex.h:31
int size() const
Get the number of attached blocks.
Definition Vertex.h:110
Pair< int > const & inPropagatorId(int i) const
Get the block and direction of an incoming propagator.
Definition Vertex.h:114
Pair< int > const & outPropagatorId(int i) const
Get the block and direction of an outgoing propagator.
Definition Vertex.h:118
Dynamically allocatable contiguous array template.
Definition DArray.h:32
void allocate(int capacity)
Allocate the underlying C array.
Definition DArray.h:199
Dynamically allocated Matrix.
Definition DMatrix.h:25
void allocate(int capacity1, int capacity2)
Allocate memory for a matrix.
Definition DMatrix.h:170
A fixed size (static) contiguous array template.
Definition FArray.h:47
An array of exactly 2 objects.
Definition Pair.h:24
ScalarParam< Type > & read(std::istream &in, const char *label, Type &value)
Add and read a new required ScalarParam < Type > object.
void setClassName(const char *className)
Set class name string.
ScalarParam< Type > & readOptional(std::istream &in, const char *label, Type &value)
Add and read a new optional ScalarParam < Type > object.
#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
bool isThread()
Is the thread model in use ?
bool isBead()
Is the bead model in use ?
PSCF package top-level namespace.
Definition param_pc.dox:1
Struct containing an enumeration of polymer structure types.
Definition PolymerType.h:27