PSCF v1.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 <prdc/cuda/VecOp.h>
17#include <pscf/math/IntVec.h>
18#include <util/random/Random.h>
19
20namespace Pscf {
21namespace Rpg {
22
23 using namespace Util;
24 using namespace Pscf::Prdc;
25 using namespace Pscf::Prdc::Cuda;
26
27 /*
28 * Constructor.
29 */
30 template <int D>
32 : BdStep<D>(simulator),
33 wp_(),
34 wf_(),
35 dci_(),
36 eta_(),
37 dwc_(),
38 dwp_(),
39 mobility_(0.0)
40 {}
41
42 /*
43 * Destructor, empty default implementation.
44 */
45 template <int D>
48
49 /*
50 * ReadParameters, empty default implementation.
51 */
52 template <int D>
53 void PredCorrBdStep<D>::readParameters(std::istream &in)
54 {
55 read(in, "mobility", mobility_);
56
57 // Allocate memory for private containers
58 int nMonomer = system().mixture().nMonomer();
59 IntVec<D> meshDimensions = system().domain().mesh().dimensions();
60 wp_.allocate(nMonomer);
61 wf_.allocate(nMonomer);
62 for (int i=0; i < nMonomer; ++i) {
63 wp_[i].allocate(meshDimensions);
64 wf_[i].allocate(meshDimensions);
65 }
66 dci_.allocate(nMonomer-1);
67 eta_.allocate(nMonomer-1);
68 for (int i=0; i < nMonomer - 1; ++i) {
69 dci_[i].allocate(meshDimensions);
70 eta_[i].allocate(meshDimensions);
71 }
72 dwc_.allocate(meshDimensions);
73 dwp_.allocate(meshDimensions);
74 }
75
76 template <int D>
78 {
79 // Check array capacities
80 int meshSize = system().domain().mesh().size();
81 int nMonomer = system().mixture().nMonomer();
82 UTIL_CHECK(wp_.capacity() == nMonomer);
83 UTIL_CHECK(wf_.capacity() == nMonomer);
84 for (int i=0; i < nMonomer; ++i) {
85 UTIL_CHECK(wp_[i].capacity() == meshSize);
86 UTIL_CHECK(wf_[i].capacity() == meshSize);
87 }
88 UTIL_CHECK(dci_.capacity() == nMonomer-1);
89 UTIL_CHECK(eta_.capacity() == nMonomer-1);
90 for (int i=0; i < nMonomer - 1; ++i) {
91 UTIL_CHECK(dci_[i].capacity() == meshSize);
92 UTIL_CHECK(eta_[i].capacity() == meshSize);
93 }
94 UTIL_CHECK(dwc_.capacity() == meshSize);
95 UTIL_CHECK(dwp_.capacity() == meshSize);
96 }
97
98 template <int D>
100 {
101 // Array sizes and indices
102 const int nMonomer = system().mixture().nMonomer();
103 const int meshSize = system().domain().mesh().size();
104 int i, j;
105
106 // Save current state
107 simulator().saveState();
108
109 // Copy current W fields from parent system
110 for (i = 0; i < nMonomer; ++i) {
111 // wp_[i] and wf_[i] set equal to w-field i
112 VecOp::eqVPair(wp_[i], wf_[i], system().w().rgrid(i));
113 }
114
115 // Store initial value of pressure field
116 VecOp::eqV(dwp_, simulator().wc(nMonomer-1));
117
118 // Constants for dynamics
119 const double vSystem = system().domain().unitCell().volume();
120 double a = -1.0*mobility_;
121 double b = sqrt(2.0*mobility_*double(meshSize)/vSystem);
122
123 // Constants for normal distribution
124 double stddev = 1.0;
125 double mean = 0;
126
127 // Construct all random displacement components
128 for (j = 0; j < nMonomer - 1; ++j) {
129 cudaRandom().normal(eta_[j], stddev, mean);
130 VecOp::mulEqS(eta_[j], b);
131 }
132
133 // Compute predicted state wp_, and store initial force dci_
134
135 // Loop over eigenvectors of projected chi matrix
136 double evec;
137 for (j = 0; j < nMonomer - 1; ++j) {
138 RField<D> const & dc = simulator().dc(j);
139 RField<D> const & eta = eta_[j];
140 RField<D> & dci = dci_[j];
141
142 // dwc_[k] = a*dc[k] + eta[k];
143 VecOp::addVcVc(dwc_, dc, a, eta, 1.0);
144
145 // dci[k] = dc[k];
146 VecOp::eqV(dci, dc);
147
148 // Loop over monomer types
149 for (i = 0; i < nMonomer; ++i) {
150 evec = simulator().chiEvecs(j,i);
151 VecOp::addEqVc(wp_[i], dwc_, evec);
152 }
153 }
154
155 // Set modified fields at predicted state wp_
156 system().w().setRGrid(wp_);
157
158 // Enforce incompressibility (also solves MDE repeatedly)
159 bool isConverged = false;
160 int compress = simulator().compressor().compress();
161 if (compress != 0){
162 simulator().restoreState();
163 } else {
164 UTIL_CHECK(system().c().hasData());
165
166 // Compute components and derivatives at wp_
167 simulator().clearData();
168 simulator().computeWc();
169 simulator().computeCc();
170 simulator().computeDc();
171
172 // Compute change in pressure field
173 RField<D> const & wp = simulator().wc(nMonomer-1);
174
175 // dwp_[k] = wp[k] - dwp_[k]
176 VecOp::subVV(dwp_, wp, dwp_);
177
178 // Adjust pressure field
179 for (i = 0; i < nMonomer; ++i) {
180 VecOp::addEqV(wf_[i], dwp_);
181 }
182
183 // Full step (corrector)
184 double ha = 0.5*a;
185 for (j = 0; j < nMonomer - 1; ++j) {
186 RField<D> const & dcp = simulator().dc(j);
187 RField<D> const & dci = dci_[j];
188 RField<D> const & eta = eta_[j];
189
190 // dwc_[k] = ha*( dci[k] + dcp[k]) + eta[k];
191 VecOp::addVcVcVc(dwc_, dci, ha, dcp, ha, eta, 1.0);
192
193 for (i = 0; i < nMonomer; ++i) {
194 evec = simulator().chiEvecs(j,i);
195 VecOp::addEqVc(wf_[i], dwc_, evec);
196 }
197 }
198
199 // Set fields at final point
200 system().w().setRGrid(wf_);
201
202 // Enforce incompressibility for final point
203 int compress2 = simulator().compressor().compress();
204 if (compress2 != 0){
205 simulator().restoreState();
206 } else {
207 isConverged = true;
208 UTIL_CHECK(system().c().hasData());
209
210 // Compute components and derivatives at final point
211 simulator().clearState();
212 simulator().clearData();
213 simulator().computeWc();
214 simulator().computeCc();
215 simulator().computeDc();
216
217 }
218 }
219
220 return isConverged;
221 }
222
223}
224}
225#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 addEqV(DeviceArray< cudaReal > &a, DeviceArray< cudaReal > const &b, const int beginIdA, const int beginIdB, const int n)
Vector addition in-place, a[i] += b[i], kernel wrapper (cudaReal).
Definition VecOp.cu:1672
void eqV(DeviceArray< cudaReal > &a, DeviceArray< cudaReal > const &b, const int beginIdA, const int beginIdB, const int n)
Vector assignment, a[i] = b[i], kernel wrapper (cudaReal).
Definition VecOp.cu:1039
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 subVV(DeviceArray< cudaReal > &a, DeviceArray< cudaReal > const &b, DeviceArray< cudaReal > const &c, const int beginIdA, const int beginIdB, const int beginIdC, const int n)
Vector subtraction, a[i] = b[i] - c[i], kernel wrapper (cudaReal).
Definition VecOp.cu:1247
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 mulEqS(DeviceArray< cudaReal > &a, const cudaReal b, const int beginIdA, const int n)
Vector multiplication in-place, a[i] *= b, kernel wrapper (cudaReal).
Definition VecOp.cu:1918
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 Reduce.cpp:14
Periodic fields and crystallography.
Definition CField.cpp:11
SCFT and PS-FTS with real periodic fields (GPU)
PSCF package top-level namespace.
Definition param_pc.dox:1