PSCF v1.3
rpc/solvers/Block.tpp
1#ifndef RPC_BLOCK_TPP
2#define RPC_BLOCK_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 "Block.h"
12
13#include <prdc/cpu/WaveList.h>
14#include <prdc/cpu/FFT.h>
15#include <prdc/crystal/UnitCell.h>
16#include <prdc/crystal/shiftToMinimum.h>
17
18#include <pscf/chem/Edge.h>
19#include <pscf/mesh/Mesh.h>
20#include <pscf/mesh/MeshIterator.h>
21#include <pscf/math/IntVec.h>
22
23#include <util/containers/DMatrix.h>
24#include <util/containers/DArray.h>
25#include <util/containers/FSArray.h>
26
27namespace Pscf {
28namespace Rpc {
29
30 // Namespaces from which names may be used without qualification
31 using namespace Util;
32 using namespace Pscf::Prdc;
33 using namespace Pscf::Prdc::Cpu;
34
35 /*
36 * Constructor.
37 */
38 template <int D>
40 : meshPtr_(nullptr),
41 fftPtr_(nullptr),
42 unitCellPtr_(nullptr),
43 waveListPtr_(nullptr),
44 kMeshDimensions_(-1),
45 kSize_(-1),
46 ds_(-1.0),
47 dsTarget_(-1.0),
48 ns_(-1),
49 isAllocated_(false),
50 hasExpKsq_(false)
51 {
52 propagator(0).setBlock(*this);
53 propagator(1).setBlock(*this);
54 }
55
56 /*
57 * Destructor.
58 */
59 template <int D>
62
63 /*
64 * Store addresses of mesh, FFT and unit cell.
65 */
66 template <int D>
67 void Block<D>::associate(Mesh<D> const & mesh,
68 FFT<D> const& fft,
69 UnitCell<D> const& cell,
70 WaveList<D>& wavelist)
71 {
72 // Preconditions
73 UTIL_CHECK(!isAllocated_);
74
75 // Set pointers to mesh and fft
76 meshPtr_ = &mesh;
77 fftPtr_ = &fft;
78 unitCellPtr_ = &cell;
79 waveListPtr_ = &wavelist;
80
81 hasExpKsq_ = false;
82 }
83
84 /*
85 * Compute number of contour steps and allocate all memory.
86 */
87 template <int D>
89 {
90 UTIL_CHECK(ds > 0.0);
91 UTIL_CHECK(meshPtr_);
92 UTIL_CHECK(fftPtr_);
93 UTIL_CHECK(unitCellPtr_);
94 UTIL_CHECK(waveListPtr_);
95 UTIL_CHECK(mesh().size() > 1);
96 UTIL_CHECK(mesh().dimensions() == fft().meshDimensions());
97 UTIL_CHECK(!isAllocated_);
98
99 // Compute DFT grid dimensions kMeshDimensions_ and size kSize_
100 FFT<D>::computeKMesh(mesh().dimensions(), kMeshDimensions_, kSize_);
101
102 // Allocate work arrays for MDE solution
103 expKsq_.allocate(kMeshDimensions_);
104 expKsq2_.allocate(kMeshDimensions_);
105 expW_.allocate(mesh().dimensions());
106 qr_.allocate(mesh().dimensions());
107 qk_.allocate(mesh().dimensions());
108 qr2_.allocate(mesh().dimensions());
109 qk2_.allocate(mesh().dimensions());
111 expW2_.allocate(mesh().dimensions());
112 } else {
113 expWInv_.allocate(mesh().dimensions());
114 }
115
116 // Allocate block concentration field
117 cField().allocate(mesh().dimensions());
118
119 dsTarget_ = ds;
120
121 // Compute ns_
123 // Set contour length discretization for this block
124 UTIL_CHECK(length() > 0.0);
125 int tempNs;
126 tempNs = floor(length() / (2.0 *ds) + 0.5);
127 if (tempNs == 0) {
128 tempNs = 1; // ensure ns_ >= 3
129 }
130 ns_ = 2*tempNs + 1;
131 ds_ = length()/double(ns_-1);
132 } else
133 if (PolymerModel::isBead()) {
134 ds_ = 1.0;
135 ns_ = nBead() + 2;
136 }
137
138 // Allocate memory for solutions to MDE (requires ns_)
139 propagator(0).allocate(ns_, mesh());
140 propagator(1).allocate(ns_, mesh());
141
142 isAllocated_ = true;
143 hasExpKsq_ = false;
144 }
145
146 /*
147 * Set or reset the the block length.
148 */
149 template <int D>
150 void Block<D>::setLength(double newLength)
151 {
153 Edge::setLength(newLength);
154
155 if (isAllocated_) {
156
157 int oldNs = ns_;
158
159 // Reset contour length discretization
160 UTIL_CHECK(dsTarget_ > 0);
161 int tempNs;
162 tempNs = floor( length()/(2.0 *dsTarget_) + 0.5 );
163 if (tempNs == 0) {
164 tempNs = 1;
165 }
166 ns_ = 2*tempNs + 1;
167 ds_ = length()/double(ns_-1);
168
169 // Reallocate propagators if ns_ has changed
170 if (oldNs != ns_) {
171 propagator(0).reallocate(ns_);
172 propagator(1).reallocate(ns_);
173 }
174
175 }
176 hasExpKsq_ = false;
177 }
178
179 /*
180 * Set or reset the the block length.
181 */
182 template <int D>
184 {
186 hasExpKsq_ = false;
187 }
188
189 /*
190 * Clear all internal data that depends on the unit cell parameters.
191 */
192 template <int D>
194 {
195 hasExpKsq_ = false;
196 stress_.clear();
197 }
198
199 /*
200 * Compute all elements of expKsq_ and expKsq2_ arrays
201 */
202 template <int D>
203 void Block<D>::computeExpKsq()
204 {
205 UTIL_CHECK(isAllocated_);
206 UTIL_CHECK(waveListPtr_);
207
208 bool isThread = PolymerModel::isThread();
209 double bSqFactor = -1.0*kuhn()*kuhn() / 6.0;
210 if (isThread) {
211 bSqFactor *= ds_;
212 }
213
214 // Calculate KSq if necessary
215 if (!waveListPtr_->hasKSq()) {
216 waveListPtr_->computeKSq();
217 }
218 RField<D> const & kSq = waveListPtr_->kSq();
219
220 MeshIterator<D> iter;
221 iter.setDimensions(kMeshDimensions_);
222 double arg;
223 int i;
224 for (iter.begin(); !iter.atEnd(); ++iter) {
225 i = iter.rank();
226 arg = kSq[i]*bSqFactor;
227 expKsq_[i] = exp(arg);
228 expKsq2_[i] = exp(0.5*arg);
229 }
230
231 hasExpKsq_ = true;
232 }
233
234 /*
235 * Setup the the step algorithm for a specific field configuration.
236 */
237 template <int D>
238 void
240 {
241 // Preconditions
242 int nx = mesh().size();
243 UTIL_CHECK(nx > 0);
244 UTIL_CHECK(isAllocated_);
245
246 // Compute expW arrays
247 double arg;
249 double c = -0.5*ds_;
250 for (int i = 0; i < nx; ++i) {
251 arg = c*w[i];
252 expW_[i] = exp(arg);
253 expW2_[i] = exp(0.5*arg);
254 }
255 } else
256 if (PolymerModel::isBead()) {
257 for (int i = 0; i < nx; ++i) {
258 arg = -w[i];
259 expW_[i] = exp(arg);
260 expWInv_[i] = 1.0/expW_[i];
261 }
262 }
263
264 // Compute expKsq arrays if necessary
265 if (!hasExpKsq_) {
266 computeExpKsq();
267 }
268
269 }
270
271 /*
272 * Propagate solution by one step for the thread model.
273 */
274 template <int D>
275 void Block<D>::stepThread(RField<D> const & q, RField<D>& qout) const
276 {
278
279 // Preconditions on mesh and fft
280 int nx = mesh().size();
281 UTIL_CHECK(nx > 0);
282 UTIL_CHECK(fft().isSetup());
283 UTIL_CHECK(mesh().dimensions() == fft().meshDimensions());
284
285 // Internal preconditions
286 UTIL_CHECK(isAllocated_);
287 int nk = qk_.capacity();
288 UTIL_CHECK(nk > 0);
289 UTIL_CHECK(qr_.capacity() == nx);
290 UTIL_CHECK(expW_.capacity() == nx);
291 UTIL_CHECK(expW2_.capacity() == nx);
292 UTIL_CHECK(expKsq_.capacity() == nk);
293 UTIL_CHECK(expKsq2_.capacity() == nk);
294 UTIL_CHECK(hasExpKsq_);
295
296 // Preconditions on parameters
298 UTIL_CHECK(q.capacity() == nx);
299 UTIL_CHECK(qout.isAllocated());
300 UTIL_CHECK(qout.capacity() == nx);
301
302 // Apply pseudo-spectral algorithm
303
304 // Step by ds/2 for qr_, step by ds/4 for qr2_
305 int i;
306 for (i = 0; i < nx; ++i) {
307 qr_[i] = q[i]*expW_[i];
308 qr2_[i] = q[i]*expW2_[i];
309 }
310 fft().forwardTransform(qr_, qk_);
311 fft().forwardTransform(qr2_, qk2_);
312 for (i = 0; i < nk; ++i) {
313 qk_[i][0] *= expKsq_[i];
314 qk_[i][1] *= expKsq_[i];
315 qk2_[i][0] *= expKsq2_[i];
316 qk2_[i][1] *= expKsq2_[i];
317 }
318 fft().inverseTransformUnsafe(qk_, qr_); // overwrites qk_
319 fft().inverseTransformUnsafe(qk2_, qr2_); // overwrites qk2_
320 for (i = 0; i < nx; ++i) {
321 qr_[i] = qr_[i]*expW_[i];
322 qr2_[i] = qr2_[i]*expW_[i];
323 }
324
325 // Above, multiplying qr2_ by expW_, rather than by expW2_, combines
326 // required multiplications by expW2_ at the end of first half-step
327 // and at the beginning of the second.
328
329 // Finish second half-step of ds/2 for qr2_
330 fft().forwardTransform(qr2_, qk2_);
331 for (i = 0; i < nk; ++i) {
332 qk2_[i][0] *= expKsq2_[i];
333 qk2_[i][1] *= expKsq2_[i];
334 }
335 fft().inverseTransformUnsafe(qk2_, qr2_); // overwrites qk2_
336 for (i = 0; i < nx; ++i) {
337 qr2_[i] = qr2_[i]*expW2_[i];
338 }
339
340 // Richardson extrapolation
341 for (i = 0; i < nx; ++i) {
342 qout[i] = (4.0*qr2_[i] - qr_[i])/3.0;
343 }
344
345 }
346
347 /*
348 * Apply one step of MDE solution for the bead model.
349 */
350 template <int D>
351 void Block<D>::stepBead(RField<D> const & q, RField<D>& qout) const
352 {
354 stepBondBead(q, qout);
355 stepFieldBead(qout);
356 }
357
358 /*
359 * Apply the bond operator for the bead model.
360 */
361 template <int D>
362 void Block<D>::stepBondBead(RField<D> const & q, RField<D>& qout) const
363 {
364 // Prereconditions
365 UTIL_CHECK(isAllocated_);
366 UTIL_CHECK(hasExpKsq_);
367 int nx = mesh().size();
368 UTIL_CHECK(nx > 0);
369 UTIL_CHECK(q.capacity() == nx);
370 UTIL_CHECK(qout.capacity() == nx);
371 int nk = qk_.capacity();
372 UTIL_CHECK(nk > 0);
373 UTIL_CHECK(expKsq_.capacity() == nk);
374
375 // Apply bond operator
376 fft().forwardTransform(q, qk_);
377 for (int i = 0; i < nk; ++i) {
378 qk_[i][0] *= expKsq_[i];
379 qk_[i][1] *= expKsq_[i];
380 }
381 fft().inverseTransformUnsafe(qk_, qout); // destroys qk_
382 }
383
384 /*
385 * Apply the half-bond operator for the bead model.
386 */
387 template <int D>
388 void Block<D>::stepHalfBondBead(RField<D> const & q, RField<D>& qout) const
389 {
390 // Preconditions
391 UTIL_CHECK(isAllocated_);
392 UTIL_CHECK(hasExpKsq_);
393 int nx = mesh().size();
394 UTIL_CHECK(nx > 0);
395 UTIL_CHECK(q.capacity() == nx);
396 UTIL_CHECK(qout.capacity() == nx);
397 int nk = qk_.capacity();
398 UTIL_CHECK(nk > 0);
399 UTIL_CHECK(expKsq_.capacity() == nk);
400
401 // Apply bond operator
402 fft().forwardTransform(q, qk_);
403 for (int i = 0; i < nk; ++i) {
404 qk_[i][0] *= expKsq2_[i];
405 qk_[i][1] *= expKsq2_[i];
406 }
407 fft().inverseTransformUnsafe(qk_, qout); // destroys qk_
408 }
409
410 /*
411 * Apply the local field operator for the bead model.
412 */
413 template <int D>
415 {
416 // Preconditions
417 int nx = mesh().size();
418 UTIL_CHECK(nx > 0);
419 UTIL_CHECK(expW_.capacity() == nx);
420 UTIL_CHECK(q.capacity() == nx);
421
422 // Apply field operator
423 for (int i = 0; i < nx; ++i) {
424 q[i] *= expW_[i];
425 }
426 }
427
428 /*
429 * Integrate to calculate monomer concentration for this block
430 */
431 template <int D>
433 {
434 // Preconditions
435 UTIL_CHECK(isAllocated_);
436 int nx = mesh().size();
437 UTIL_CHECK(nx > 0);
438 UTIL_CHECK(ns_ > 0);
439 UTIL_CHECK(propagator(0).isSolved());
440 UTIL_CHECK(propagator(1).isSolved());
441 UTIL_CHECK(cField().capacity() == nx);
442
443 // Initialize cField to zero at all points
444 int i;
445 for (i = 0; i < nx; ++i) {
446 cField()[i] = 0.0;
447 }
448
449 // References to forward and reverse propagators
450 Propagator<D> const & p0 = propagator(0);
451 Propagator<D> const & p1 = propagator(1);
452
453 // Evaluate unnormalized integral
454
455 // Endpoint contributions
456 for (i = 0; i < nx; ++i) {
457 cField()[i] += p0.q(0)[i]*p1.q(ns_ - 1)[i];
458 cField()[i] += p0.q(ns_ -1)[i]*p1.q(0)[i];
459 }
460
461 // Odd indices
462 int j;
463 for (j = 1; j < (ns_ -1); j += 2) {
464 RField<D> const & qf = p0.q(j);
465 RField<D> const & qr = p1.q(ns_ - 1 - j);
466 for (i = 0; i < nx; ++i) {
467 //cField()[i] += p0.q(j)[i] * p1.q(ns_ - 1 - j)[i] * 4.0;
468 cField()[i] += qf[i] * qr[i] * 4.0;
469 }
470 }
471
472 // Even indices
473 for (j = 2; j < (ns_ -2); j += 2) {
474 RField<D> const & qf = p0.q(j);
475 RField<D> const & qr = p1.q(ns_ - 1 - j);
476 for (i = 0; i < nx; ++i) {
477 // cField()[i] += p0.q(j)[i] * p1.q(ns_ - 1 - j)[i] * 2.0;
478 cField()[i] += qf[i] * qr[i] * 2.0;
479 }
480 }
481
482 // Normalize the integral
483 prefactor *= ds_ / 3.0;
484 for (i = 0; i < nx; ++i) {
485 cField()[i] *= prefactor;
486 }
487
488 }
489
490 /*
491 * Calculate monomer concentration for this block, bead model.
492 */
493 template <int D>
495 {
496 // Preconditions
497 UTIL_CHECK(isAllocated_);
498 int nx = mesh().size();
499 UTIL_CHECK(nx > 0);
500 UTIL_CHECK(ns_ > 0);
501 UTIL_CHECK(propagator(0).isSolved());
502 UTIL_CHECK(propagator(1).isSolved());
503 UTIL_CHECK(cField().capacity() == nx);
504
505 // Initialize cField to zero at all points
506 int i;
507 for (i = 0; i < nx; ++i) {
508 cField()[i] = 0.0;
509 }
510
511 // References to forward and reverse propagators
512 Propagator<D> const & p0 = propagator(0);
513 Propagator<D> const & p1 = propagator(1);
514
515 // Sum over beads (j = 1, ... , ns_ -2)
516 int j;
517 for (j = 1; j < (ns_ -1); ++j) {
518 RField<D> const & qf = p0.q(j);
519 RField<D> const & qr = p1.q(ns_ - 1 - j);
520 for (i = 0; i < nx; ++i) {
521 cField()[i] += qf[i] * qr[i] * expWInv_[i];
522 }
523 }
524
525 // Normalize the integral
526 for (i = 0; i < nx; ++i) {
527 cField()[i] *= prefactor;
528 }
529 }
530
531 /*
532 * Integrate to compute stress exerted by this block (thread).
533 */
534 template <int D>
535 void Block<D>::computeStressThread(double prefactor)
536 {
537 // Preconditions
539 int nx = mesh().size();
540 UTIL_CHECK(nx > 0);
541 UTIL_CHECK(fft().isSetup());
542 UTIL_CHECK(mesh().dimensions() == fft().meshDimensions());
543 UTIL_CHECK(ns_ > 0);
544 UTIL_CHECK(ds_ > 0);
545
546 // References to forward and reverse propagators
547 Propagator<D> const & p0 = propagator(0);
548 Propagator<D> const & p1 = propagator(1);
551
552 // If necessary, update derivatives of |k|^2 in WaveList
553 UTIL_CHECK(waveListPtr_);
554 if (!waveListPtr_->hasdKSq()) {
555 waveListPtr_->computedKSq();
556 }
557
558 // Initialize dQ work array to zero
560 const int nParam = unitCell().nParameter();
561 for (int i = 0; i < nParam; ++i) {
562 dQ.append(0.0);
563 }
564
565 // Work variables
566 const double bSq = kuhn()*kuhn()/6.0;
567 double dels, prod, increment;
568 int m, n;
569
570 // Evaluate unnormalized integral over contour
571 for (int j = 0; j < ns_ ; ++j) {
572
573 qr_ = p0.q(j);
574 fft().forwardTransform(qr_, qk_);
575
576 qr2_ = p1.q(ns_ - 1 - j);
577 fft().forwardTransform(qr2_, qk2_);
578
579 // Compute prefactor dels for Simpson's rule
580 dels = ds_ / 3.0;
581 if (j != 0 && j != ns_ - 1) {
582 if (j % 2 == 0) {
583 dels *= 2.0;
584 } else {
585 dels *= 4.0;
586 }
587 }
588
589 // Loop over unit cell parameters
590 for (n = 0; n < nParam ; ++n) {
591 RField<D> dKSq = waveListPtr_->dKSq(n);
592 increment = 0.0;
593 // Loop over wavevectors
594 for (m = 0; m < kSize_ ; ++m) {
595 prod = (qk2_[m][0] * qk_[m][0]) + (qk2_[m][1] * qk_[m][1]);
596 prod *= dKSq[m];
597 increment += prod;
598 }
599 increment *= bSq * dels;
600 dQ[n] = dQ[n] - increment;
601 }
602
603 }
604
605 // Compute stress_ from dQ
606 stress_.clear();
607 for (int n = 0; n < nParam; ++n) {
608 stress_.append(-1.0*prefactor*dQ[n]);
609 }
610
611 }
612
613 /*
614 * Compute contribution of this block to stress for bead model.
615 */
616 template <int D>
617 void Block<D>::computeStressBead(double prefactor)
618 {
619 // Preconditions
621 int nx = mesh().size();
622 UTIL_CHECK(nx > 0);
623 UTIL_CHECK(fft().isSetup());
624 UTIL_CHECK(mesh().dimensions() == fft().meshDimensions());
625 UTIL_CHECK(ns_ > 0);
626
627 // References to forward to reverse propagators
628 Propagator<D> const & p0 = propagator(0);
629 Propagator<D> const & p1 = propagator(1);
630 UTIL_CHECK(p0.isSolved());
631 UTIL_CHECK(p1.isSolved());
632
633
634 // If necessary, update derivatives of |k|^2 in WaveList
635 UTIL_CHECK(waveListPtr_);
636 if (!waveListPtr_->hasdKSq()) {
637 waveListPtr_->computedKSq();
638 }
639
640 // Initialize dQ work array to zero
642 int nParam = unitCell().nParameter();
643 for (int i = 0; i < nParam; ++i) {
644 dQ.append(0.0);
645 }
646
647 const double bSq = kuhn()*kuhn()/6.0;
648 double increment, prod;
649
650 // Half-bond from head junction to first bead, if not a chain end
651 if (!p0.isHeadEnd()) {
652
653 // Bead j = 0, forward propagator
654 qr_ = p0.q(0);
655 fft().forwardTransform(qr_, qk_);
656
657 // Next bead (ns_ - 2), reverse propagator
658 qr2_ = p1.q(ns_ - 2);
659 fft().forwardTransform(qr2_, qk2_);
660
661 // Loop over unit cell parameters
662 for (int n = 0; n < nParam ; ++n) {
663 RField<D> dKSq = waveListPtr_->dKSq(n);
664 increment = 0.0;
665 // Loop over wavevectors
666 for (int m = 0; m < kSize_ ; ++m) {
667 prod = (qk2_[m][0] * qk_[m][0]) + (qk2_[m][1] * qk_[m][1]);
668 prod *= dKSq[m]*expKsq2_[m];
669 increment += prod;
670 }
671 increment *= 0.5*bSq;
672 dQ[n] = dQ[n] - increment;
673 }
674 }
675
676 // Loop over nBead - 1 full bonds within this block
677 for (int j = 1; j < ns_ - 2 ; ++j) {
678
679 // Bead j, forward propagator
680 qr_ = p0.q(j);
681 fft().forwardTransform(qr_, qk_);
682
683 // Next bead (ns_- 2 - j), reverse propagator
684 qr2_ = p1.q(ns_ - 2 - j);
685 fft().forwardTransform(qr2_, qk2_);
686
687 // Loop over unit cell parameters
688 for (int n = 0; n < nParam ; ++n) {
689 RField<D> dKSq = waveListPtr_->dKSq(n);
690 increment = 0.0;
691 // Loop over wavevectors
692 for (int m = 0; m < kSize_ ; ++m) {
693 prod = (qk2_[m][0] * qk_[m][0]) + (qk2_[m][1] * qk_[m][1]);
694 prod *= dKSq[m]*expKsq_[m];
695 increment += prod;
696 }
697 increment *= bSq;
698 dQ[n] = dQ[n] - increment;
699 }
700
701 }
702
703 // Half-bond from last bead to tail junction, if not a chain end
704 if (!p0.isTailEnd()) {
705
706 // Second-to-last bead, j = ns_ - 2, forward propagator
707 qr_ = p0.q(ns_-2);
708 fft().forwardTransform(qr_, qk_);
709
710 // Last bead, j=0, reverse propagator
711 qr2_ = p1.q(0);
712 fft().forwardTransform(qr2_, qk2_);
713
714 // Loop over unit cell parameters
715 for (int n = 0; n < nParam ; ++n) {
716 RField<D> dKSq = waveListPtr_->dKSq(n);
717 increment = 0.0;
718 // Loop over wavevectors
719 for (int m = 0; m < kSize_ ; ++m) {
720 prod = (qk2_[m][0] * qk_[m][0]) + (qk2_[m][1] * qk_[m][1]);
721 prod *= dKSq[m]*expKsq2_[m];
722 increment += prod;
723 }
724 increment *= 0.5*bSq;
725 dQ[n] = dQ[n] - increment;
726 }
727
728 }
729
730 // Compute stress_ from dQ
731 stress_.clear();
732 for (int i = 0; i < nParam; ++i) {
733 stress_.append(-1.0*prefactor*dQ[i]);
734 }
735
736 }
737
738}
739}
740#endif
virtual void setLength(double length)
Set the length of this block (only valid for thread model).
Definition Edge.cpp:68
Iterator over points in a Mesh<D>.
int rank() const
Get the rank of current element.
void begin()
Set iterator to the first point in the mesh.
bool atEnd() const
Is this the end (i.e., one past the last point)?
void setDimensions(const IntVec< D > &dimensions)
Set the grid dimensions in all directions.
Description of a regular grid of points in a periodic domain.
Definition Mesh.h:61
Fourier transform wrapper.
Definition cpu/FFT.h:38
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
bool isAllocated() const
Return true if the FftwDArray has been allocated, false otherwise.
Definition FftwDArray.h:102
Field of real double precision values on an FFT mesh.
Definition cpu/RField.h:29
Class to calculate and store properties of wavevectors.
Base template for UnitCell<D> classes, D=1, 2 or 3.
Definition UnitCell.h:56
void stepBondBead(RField< D > const &qin, RField< D > &qout) const
Apply a bond operator for the bead model.
void setLength(double newLength)
Set or reset block length (only used in thread model).
void computeConcentrationBead(double prefactor)
Compute the concentration for this block, using the bead model.
void setKuhn(double kuhn)
Set or reset monomer statistical segment length.
void clearUnitCellData()
Clear all internal data that depends on the unit cell parameters.
void setupSolver(RField< D > const &w)
Set up the MDE solver for this block.
void allocate(double ds)
Allocate memory and set contour step size.
void computeStressBead(double prefactor)
Compute stress contribution for this block, using bead model.
int nBead() const
Get the number of beads in this block, in the bead model.
Definition Edge.h:302
Propagator< D > & propagator(int directionId)
Get a Propagator for a specified direction.
void associate(Mesh< D > const &mesh, FFT< D > const &fft, UnitCell< D > const &cell, WaveList< D > &wavelist)
Create permanent associations with related objects.
double kuhn() const
Get monomer statistical segment length.
Propagator< D >::FieldT & cField()
Get the associated monomer concentration field.
void computeConcentrationThread(double prefactor)
Compute the concentration for this block, for the thread model.
void stepThread(RField< D > const &qin, RField< D > &qout) const
Compute one step of solution of MDE for the thread model.
void computeStressThread(double prefactor)
Compute stress contribution for this block, using thread model.
double ds() const
Get contour length step size.
void stepFieldBead(RField< D > &q) const
Apply the exponential field operator for the bead model.
void stepHalfBondBead(RField< D > const &qin, RField< D > &qout) const
Apply a half-bond operator for the bead model.
void stepBead(RField< D > const &qin, RField< D > &qout) const
Compute one step of solution of MDE for the bead model.
double length() const
Get the length of this block, in the thread model.
Definition Edge.h:311
MDE solver for one direction of one block.
bool isHeadEnd() const
Is the head vertex a chain end?
bool isSolved() const
Has the modified diffusion equation been solved?
bool isTailEnd() const
Is the tail vertex a chain end?
bool isAllocated() const
Has memory been allocated for this propagator?
const FieldT & q(int i) const
Return q-field at specified step.
int capacity() const
Return allocated size.
Definition Array.h:159
A fixed capacity (static) contiguous array with a variable logical size.
Definition FSArray.h:38
void append(Data const &data)
Append data to the end of the array.
Definition FSArray.h:258
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
Definition global.h:68
bool isThread()
Is the thread model in use ?
bool isBead()
Is the bead model in use ?
Fields and FFTs for periodic boundary conditions (CPU)
Definition CField.cpp:12
Periodic fields and crystallography.
Definition CField.cpp:11
Real periodic fields, SCFT and PS-FTS (CPU).
Definition param_pc.dox:2
PSCF package top-level namespace.
Definition param_pc.dox:1