PSCF v1.3.3
rpg/fts/brownian/PredCorrBdStep.tpp
1#ifndef RPG_PRED_CORR_BD_STEP_TPP
2#define RPG_PRED_CORR_BD_STEP_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 "PredCorrBdStep.h"
12
13#include <rpg/fts/brownian/BdSimulator.h>
14#include <rpg/fts/compressor/Compressor.h>
15#include <rpg/system/System.h>
16#include <rpg/solvers/Mixture.h>
17#include <rpg/field/Domain.h>
18#include <prdc/cuda/VecOp.h>
19#include <pscf/math/IntVec.h>
20#include <util/random/Random.h>
21
22namespace Pscf {
23namespace Rpg {
24
25 using namespace Util;
26 using namespace Pscf::Prdc;
27 using namespace Pscf::Prdc::Cuda;
28
29 /*
30 * Constructor.
31 */
32 template <int D>
34 : BdStep<D>(simulator),
35 wp_(),
36 wf_(),
37 dci_(),
38 eta_(),
39 dwc_(),
40 dwp_(),
41 mobility_(0.0)
42 {}
43
44 /*
45 * Destructor, empty default implementation.
46 */
47 template <int D>
50
51 /*
52 * ReadParameters, empty default implementation.
53 */
54 template <int D>
55 void PredCorrBdStep<D>::readParameters(std::istream &in)
56 {
57 read(in, "mobility", mobility_);
58
59 // Allocate memory for private containers
60 int nMonomer = system().mixture().nMonomer();
61 IntVec<D> meshDimensions = system().domain().mesh().dimensions();
62 wp_.allocate(nMonomer);
63 wf_.allocate(nMonomer);
64 for (int i=0; i < nMonomer; ++i) {
65 wp_[i].allocate(meshDimensions);
66 wf_[i].allocate(meshDimensions);
67 }
68 dci_.allocate(nMonomer-1);
69 eta_.allocate(nMonomer-1);
70 for (int i=0; i < nMonomer - 1; ++i) {
71 dci_[i].allocate(meshDimensions);
72 eta_[i].allocate(meshDimensions);
73 }
74 dwc_.allocate(meshDimensions);
75 dwp_.allocate(meshDimensions);
76 }
77
78 template <int D>
80 {
81 // Check array capacities
82 int meshSize = system().domain().mesh().size();
83 int nMonomer = system().mixture().nMonomer();
84 UTIL_CHECK(wp_.capacity() == nMonomer);
85 UTIL_CHECK(wf_.capacity() == nMonomer);
86 for (int i=0; i < nMonomer; ++i) {
87 UTIL_CHECK(wp_[i].capacity() == meshSize);
88 UTIL_CHECK(wf_[i].capacity() == meshSize);
89 }
90 UTIL_CHECK(dci_.capacity() == nMonomer-1);
91 UTIL_CHECK(eta_.capacity() == nMonomer-1);
92 for (int i=0; i < nMonomer - 1; ++i) {
93 UTIL_CHECK(dci_[i].capacity() == meshSize);
94 UTIL_CHECK(eta_[i].capacity() == meshSize);
95 }
96 UTIL_CHECK(dwc_.capacity() == meshSize);
97 UTIL_CHECK(dwp_.capacity() == meshSize);
98 }
99
100 template <int D>
102 {
103 // Array sizes and indices
104 const int nMonomer = system().mixture().nMonomer();
105 const int meshSize = system().domain().mesh().size();
106 int i, j;
107
108 // Save current state
109 simulator().saveState();
110
111 // Copy current W fields from parent system
112 for (i = 0; i < nMonomer; ++i) {
113 // wp_[i] and wf_[i] set equal to w-field i
114 VecOp::eqVPair(wp_[i], wf_[i], system().w().rgrid(i));
115 }
116
117 // Store initial value of pressure field
118 VecOp::eqV(dwp_, simulator().wc(nMonomer-1));
119
120 // Constants for dynamics
121 const double vSystem = system().domain().unitCell().volume();
122 double a = -1.0*mobility_;
123 double b = sqrt(2.0*mobility_*double(meshSize)/vSystem);
124
125 // Constants for normal distribution
126 double stddev = 1.0;
127 double mean = 0;
128
129 // Construct all random displacement components
130 for (j = 0; j < nMonomer - 1; ++j) {
131 cudaRandom().normal(eta_[j], stddev, mean);
132 VecOp::mulEqS(eta_[j], b);
133 }
134
135 // Compute predicted state wp_, and store initial force dci_
136
137 // Loop over eigenvectors of projected chi matrix
138 double evec;
139 for (j = 0; j < nMonomer - 1; ++j) {
140 RField<D> const & dc = simulator().dc(j);
141 RField<D> const & eta = eta_[j];
142 RField<D> & dci = dci_[j];
143
144 // dwc_[k] = a*dc[k] + eta[k];
145 VecOp::addVcVc(dwc_, dc, a, eta, 1.0);
146
147 // dci[k] = dc[k];
148 VecOp::eqV(dci, dc);
149
150 // Loop over monomer types
151 for (i = 0; i < nMonomer; ++i) {
152 evec = simulator().chiEvecs(j,i);
153 VecOp::addEqVc(wp_[i], dwc_, evec);
154 }
155 }
156
157 // Set modified fields at predicted state wp_
158 system().w().setRGrid(wp_);
159
160 // Enforce incompressibility (also solves MDE repeatedly)
161 bool isConverged = false;
162 int compress = simulator().compressor().compress();
163 if (compress != 0){
164 simulator().restoreState();
165 } else {
166 UTIL_CHECK(system().c().hasData());
167
168 // Compute components and derivatives at wp_
169 simulator().clearData();
170 simulator().computeWc();
171 simulator().computeCc();
172 simulator().computeDc();
173
174 // Compute change in pressure field
175 RField<D> const & wp = simulator().wc(nMonomer-1);
176
177 // dwp_[k] = wp[k] - dwp_[k]
178 VecOp::subVV(dwp_, wp, dwp_);
179
180 // Adjust pressure field
181 for (i = 0; i < nMonomer; ++i) {
182 VecOp::addEqV(wf_[i], dwp_);
183 }
184
185 // Full step (corrector)
186 double ha = 0.5*a;
187 for (j = 0; j < nMonomer - 1; ++j) {
188 RField<D> const & dcp = simulator().dc(j);
189 RField<D> const & dci = dci_[j];
190 RField<D> const & eta = eta_[j];
191
192 // dwc_[k] = ha*( dci[k] + dcp[k]) + eta[k];
193 VecOp::addVcVcVc(dwc_, dci, ha, dcp, ha, eta, 1.0);
194
195 for (i = 0; i < nMonomer; ++i) {
196 evec = simulator().chiEvecs(j,i);
197 VecOp::addEqVc(wf_[i], dwc_, evec);
198 }
199 }
200
201 // Set fields at final point
202 system().w().setRGrid(wf_);
203
204 // Enforce incompressibility for final point
205 int compress2 = simulator().compressor().compress();
206 if (compress2 != 0){
207 simulator().restoreState();
208 } else {
209 isConverged = true;
210 UTIL_CHECK(system().c().hasData());
211
212 // Compute components and derivatives at final point
213 simulator().clearState();
214 simulator().clearData();
215 simulator().computeWc();
216 simulator().computeCc();
217 simulator().computeDc();
218
219 }
220 }
221
222 return isConverged;
223 }
224
225}
226}
227#endif
An IntVec<D, T> is a D-component vector of elements of integer type T.
Definition IntVec.h:27
Field of real double precision values on an FFT mesh.
Definition cpu/RField.h:29
Brownian dynamics simulator.
BdStep(BdSimulator< D > &simulator)
Constructor.
System< D > & system()
Get parent System object.
CudaRandom & cudaRandom()
Get Random number generator of parent System.
BdSimulator< D > & simulator()
Get parent BdSimulator object.
virtual void setup()
Setup before simulation.
virtual bool step()
Take a single Brownian dynamics step.
ScalarParam< Type > & read(std::istream &in, const char *label, Type &value)
Add and read a new required ScalarParam < Type > object.
PredCorrBdStep(BdSimulator< D > &simulator)
Constructor.
virtual void readParameters(std::istream &in)
Read required parameters from file.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
Definition global.h:68
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 mulEqS(Array< double > &a, double b)
Vector multiplication in-place, a[i] *= b.
Definition VecOp.cpp:267
void addVcVcVc(DeviceArray< cudaReal > &a, DeviceArray< cudaReal > const &b, cudaReal const c, DeviceArray< cudaReal > const &d, cudaReal const e, DeviceArray< cudaReal > const &f, cudaReal const g)
3-vec addition w coeff, a[i] = (b[i]*c) + (d[i]*e) + (f[i]*g), kernel wrapper.
Definition VecOpMisc.cu:322
void addVcVc(DeviceArray< cudaReal > &a, DeviceArray< cudaReal > const &b, cudaReal const c, DeviceArray< cudaReal > const &d, cudaReal const e)
Vector addition w/ coefficient, a[i] = (b[i]*c) + (d[i]*e), kernel wrapper.
Definition VecOpMisc.cu:304
void eqVPair(DeviceArray< cudaReal > &a1, DeviceArray< cudaReal > &a2, DeviceArray< cudaReal > const &s)
Vector assignment in pairs, ax[i] = b[i], kernel wrapper.
Definition VecOpMisc.cu:409
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
Fields, FFTs, and utilities for periodic boundary conditions (CUDA)
Definition CField.cu:12
Periodic fields and crystallography.
Definition CField.cpp:11
SCFT and PS-FTS with real periodic fields (GPU)
PSCF package top-level namespace.