-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathCrankNicholsonSolver1.H
234 lines (209 loc) · 7.93 KB
/
CrankNicholsonSolver1.H
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
//////////////////////////////////////////////////////////////////////
// Copyright 2014-2016 Jeffrey Comer
//
// This file is part of DiffusionFusion.
//
// DiffusionFusion is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
//
// DiffusionFusion is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along with DiffusionFusion. If not, see http://www.gnu.org/licenses/.
///////////////////////////////////////////////////////////////////////
// Interface for a Smoluchowski equation solver
// Author: Jeff Comer <jeffcomer at gmail>
#ifndef CRANKNICHOLSONSOLVER1_H
#define CRANKNICHOLSONSOLVER1_H
#include "Piecewise1d.H"
#include <omp.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
class CrankNicholsonSolver1 {
private:
int n;
double dx, x0, x1;
bool periodic;
double timestep;
double beta;
public:
CrankNicholsonSolver1(const Piecewise1d* refField, double timestep0, double kT) {
n = refField->length();
dx = refField->getDr();
x0 = refField->getR0();
x1 = refField->getR1();
periodic = refField->getPeriodic();
timestep = timestep0;
beta = 1.0/kT;
}
int length() const { return n; }
void init(Piecewise1d* prob, int node) const {
// Zero the initial distribution.
prob->zero();
// Add the delta function at the initial point.
prob->set(node, 1.0/dx);
}
double gaussianReflect(double x, double mu, double sigma) const {
double norm = 1.0/(sqrt(2*M_PI)*sigma);
double sigmaSq = sigma*sigma;
return norm*(exp(-0.5*(x-mu)*(x-mu)/sigmaSq) + exp(-0.5*(x-2.0*x0+mu)*(x-2.0*x0+mu)/sigmaSq) + exp(-0.5*(x-2.0*x1+mu)*(x-2.0*x1+mu)/sigmaSq));
}
// void init(double* prob, int node) const {
// init(prob,node,dx);
// }
void init(Piecewise1d* prob, double x, double width) const {
for (int i = 0; i < n; i++)
prob->set(i, gaussianReflect(prob->nodePos(i), x, width));
}
void solve(Piecewise1d* prob, int steps, const Piecewise1d* diffuse, const Piecewise1d* force, const Piecewise1d* bias) const {
// Coefficients contain the diffusion and drift effects.
double ca[n];
double cb[n];
double r[n];
for (int i = 0; i < n; i++) {
// Get the diffusivity and force.
double dif, gradDif;
double frc, gradFrc;
double fb, gradFb;
double x = x0 + i*dx;
diffuse->computeValGrad(x, dif, gradDif);
force->computeValGrad(x, frc, gradFrc);
// dif = diffuse->get(i);
// gradDif = diffuse->getGrad(i);
// frc = force->get(i);
// gradFrc = force->getGrad(i);
if (bias != NULL) {
bias->computeValGrad(x, fb, gradFb);
frc += fb;
gradFrc += gradFb;
}
/////////////////////////////////////////////////////////////////
// We are solving the differential equation for p(x,t)
// p_t = D*p_xx + (D_x - beta*D*F)*p_x - beta*(D*F_x + D_x*F)*p
// Here, subscripts represent partial derivatives.
// D(x): the diffusivity
// F(x,t) = -w_x(x) + fBias(x,t): the total force
// beta: 1/(kT)
//
// We write this equation as
// p_t = ca*p + cb*p_x + D*p_xx
// Below we use r = D*dt/dx^2, consistent with
// Numerical Methods Using Matlab 4th Ed. Mathews and Fink, pp. 561-562
// Calculate the coefficients ca, cb, and r.
ca[i] = -beta*(dif*gradFrc + gradDif*frc)*timestep;
cb[i] = 0.5*(gradDif - beta*dif*frc)/dx*timestep;
r[i] = dif*timestep/(dx*dx);
}
// Boundary conditions.
// Left
double frcA = force->computeVal(x0);
if (bias != NULL) frcA += bias->computeVal(x0);
// A dimensionless variable related to the drift at the left boundary.
double driftA = beta*frcA*dx;
// Right
double x1 = x0 + (n-1)*dx;
double frcB = force->computeVal(x1);
if (bias != NULL) frcB += bias->computeVal(x1);
// A dimensionless variable related to the drift at the right boundary.
double driftB = beta*frcB*dx;
// Allocate linear system vectors.
// We have a tridiagonal matrix.
// See: https://www.gnu.org/software/gsl/manual/html_node/Tridiagonal-Systems.html
gsl_vector* d = gsl_vector_alloc(n); // Matrix diagonal
gsl_vector* e = gsl_vector_alloc(n-1); // Band to the right of diagonal
gsl_vector* f = gsl_vector_alloc(n-1); // Band to the left of diagonal
gsl_vector* b = gsl_vector_alloc(n); // the right side of Mx = b
gsl_vector* sol = gsl_vector_alloc(n); // the solution "x"
gsl_vector* sol0 = gsl_vector_alloc(n); // a buffer to swap
// Fill in the matrix.
// Robin boundary conditions (zero flux).
// Left: Use the O(dx^2) forward finite difference.
{
int i = 0;
double corrA = (cb[i]-r[i])/(3.0+2.0*driftA);
gsl_vector_set(d, i, 2.0 + 2.0*r[i] - ca[i] + 4.0*corrA);
gsl_vector_set(e, i, -(cb[i] + r[i] + corrA));
}
// Right: Use the O(dx^2) backward finite difference.
{
int i = n-1;
double corrB = -(cb[i]+r[i])/(3.0-2.0*driftB);
gsl_vector_set(d, i, 2.0 + 2.0*r[i] - ca[i] + 4.0*corrB);
gsl_vector_set(f, i-1, cb[i] - (r[i] + corrB));
}
// Interior nodes using a Crank-Nicholson approach.
// See: Numerical Methods Using Matlab 4th Ed. Mathews and Fink, pp. 561-562
// See: https://www.gnu.org/software/gsl/manual/html_node/Tridiagonal-Systems.html
for (int i = 1; i < n-1; i++) {
gsl_vector_set(d, i, 2.0 + 2.0*r[i] - ca[i]);
gsl_vector_set(f, i-1, cb[i] - r[i]);
gsl_vector_set(e, i, -(cb[i] + r[i]));
}
// Absorbing at left.
//gsl_vector_set(d, 0, 1.0);
//gsl_vector_set(e, 0, 0.0);
//gsl_vector_set(b, 0, 0.0);
// Absorbing at right.
//gsl_vector_set(f, n-2, 0.0);
//gsl_vector_set(d, n-1, 1.0);
//gsl_vector_set(b, n-1, 0.0);
// Intialize the solution vector.
for (int i = 0; i < n; i++) gsl_vector_set(sol, i, prob->get(i));
for (int s = 0; s < steps; s++) {
// Swap the pointers.
gsl_vector* tmp = sol0;
sol0 = sol;
sol = tmp;
// Fill in the vector b.
// Reflecting at left.
{
int i = 0;
double um = (4.0*gsl_vector_get(sol0,i) - gsl_vector_get(sol0,i+1))/(3.0+2.0*driftA);
gsl_vector_set(b, i,
(2.0 - 2.0*r[i] + ca[i])*gsl_vector_get(sol0,i)
+ (r[i] - cb[i])*um
+ (r[i] + cb[i])*gsl_vector_get(sol0,i+1));
}
// Reflecting at right.
{
int i = n-1;
double up = (4.0*gsl_vector_get(sol0,i) - gsl_vector_get(sol0,i-1))/(3.0-2.0*driftB);
gsl_vector_set(b, i,
(2.0 - 2.0*r[i] + ca[i])*gsl_vector_get(sol0,i)
+ (r[i] - cb[i])*gsl_vector_get(sol0,i-1)
+ (r[i] + cb[i])*up);
}
// Interior
for (int i = 1; i < n-1; i++) {
gsl_vector_set(b, i,
(2.0 - 2.0*r[i] + ca[i])*gsl_vector_get(sol0, i)
+ (r[i] - cb[i])*gsl_vector_get(sol0, i-1)
+ (r[i] + cb[i])*gsl_vector_get(sol0, i+1));
}
// Solve the system of linear equations.
gsl_linalg_solve_tridiag(d, e, f, b, sol);
}
// Copy the results into prob.
for (int i = 0; i < n; i++) prob->set(i, gsl_vector_get(sol, i));
// Deallocate everything.
gsl_vector_free(d);
gsl_vector_free(e);
gsl_vector_free(f);
gsl_vector_free(b);
gsl_vector_free(sol);
gsl_vector_free(sol0);
}
double conserveProb(Piecewise1d* sol) const {
double sum = sol->integrate(sol->getR0(), sol->getR1());
double scale = 1.0/sum;
sol->scale(scale);
return scale;
}
int positive(Piecewise1d* sol) const {
for (int i = 0; i < sol->length(); i++) {
if (sol->get(i) < -0.5/dx) return i;
}
return -1;
}
~CrankNicholsonSolver1() {
}
};
#endif