-
Notifications
You must be signed in to change notification settings - Fork 4
/
electrostatSolver.h
336 lines (286 loc) · 13.1 KB
/
electrostatSolver.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
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
/** \file electrostatSolver.h
\brief solver for electrostatic problem when STT is required
header containing electrostatSolver class. It uses biconjugate stabilized gradient with diagonal
preconditioner. The solver is only called once to compute voltages V for each nodes of the mesh,
when STT computation is involved.
*/
#include <iostream>
#include <map>
#include "config.h"
#include "fem.h"
#include "mesh.h"
#include "facette.h"
#include "tetra.h"
#include "spinTransferTorque.h"
/** assemble the matrix K from tet and Ke inputs */
inline void assembling_mat(Tetra::Tet const &tet, double Ke[Tetra::N][Tetra::N], std::vector<Eigen::Triplet<double>> &K)
{
for (int ie = 0; ie < Tetra::N; ie++)
{
for (int je = 0; je < Tetra::N; je++)
{
double val = Ke[ie][je];
if (val != 0)
{ K.push_back( Eigen::Triplet(tet.ind[ie], tet.ind[je], val) ); }
}
}
}
/** assemble the vector L from fac and Le inputs */
inline void assembling_vect(Facette::Fac const &fac, std::vector<double> const &Le, Eigen::Ref<Eigen::VectorXd> L)
{
for (int ie = 0; ie < Facette::N; ie++)
{ L(fac.ind[ie]) += Le[ie]; }
}
/** compute side problem (electrostatic potential on the nodes) integrales for matrix
* coefficients,input from tet ; sigma is the region conductivity */
void integrales(Tetra::Tet const &tet, double sigma, double AE[Tetra::N][Tetra::N])
{
for (int npi = 0; npi < Tetra::NPI; npi++)
{
double w = tet.weight[npi];
for (int ie = 0; ie < Tetra::N; ie++)
{
double dai_dx = tet.da(ie,Nodes::IDX_X);
double dai_dy = tet.da(ie,Nodes::IDX_Y);
double dai_dz = tet.da(ie,Nodes::IDX_Z);
for (int je = 0; je < Tetra::N; je++)
{
double daj_dx = tet.da(je,Nodes::IDX_X);
double daj_dy = tet.da(je,Nodes::IDX_Y);
double daj_dz = tet.da(je,Nodes::IDX_Z);
AE[ie][je] += sigma * (dai_dx * daj_dx + dai_dy * daj_dy + dai_dz * daj_dz) * w;
}
}
}
}
/** compute integrales for vector coefficients, input from facette */
void integrales(Facette::Fac const &fac, double pot_val, std::vector<double> &BE)
{
for (int npi = 0; npi < Facette::NPI; npi++)
for (int ie = 0; ie < Facette::N; ie++)
{
BE[ie] -= Facette::a[ie][npi] * pot_val * fac.weight[npi];
}
}
/** \class electrostatSolver
this class is containing both data and a solver to compute potential from dirichlet boundary
conditions problem for the current density flowing in the sample.
*/
class electrostatSolver
{
public:
/** constructor */
inline electrostatSolver(
Mesh::mesh const &_msh /**< [in] reference to the mesh */,
STT const &_p_stt /**< all spin transfer torque parameters */,
const double _tol /**< [in] tolerance for solvers */,
const bool v /**< [in] verbose bool */,
const int max_iter /**< [in] maximum number of iteration */,
const std::string
_fileName /**< [in] output .sol file name for electrostatic potential */)
: msh(_msh), p_stt(_p_stt), verbose(v), MAXITER(max_iter), precision(PRECISION_STT),
fileName(_fileName)
{
ksi = Nodes::sq(p_stt.lJ / p_stt.lsf);
D0 = 2.0 * p_stt.sigma / (Nodes::sq(CHARGE_ELECTRON) * p_stt.N0);
pf = Nodes::sq(p_stt.lJ) / (D0 * (1. + ksi * ksi)) * BOHRS_MUB * p_stt.beta / CHARGE_ELECTRON;
if (verbose)
{
std::cout << "Dirichlet boundary conditions..." << std::endl;
infos();
}
V.resize(_msh.getNbNodes());
bool has_converged = solve(_tol);
if (has_converged)
{
if (p_stt.V_file)
{
if (verbose)
{
std::cout << "writing electrostatic potential solutions to file " << fileName
<< std::endl;
}
bool iznogood = msh.savesol(precision, fileName, "## columns: index\tV\n", V);
if (verbose && iznogood)
{
std::cout << "file " << fileName << " status : " << iznogood << std::endl;
}
}
std::for_each(msh.tet.begin(), msh.tet.end(),
[this](Tetra::Tet const &tet)
{
Eigen::Matrix<double,Nodes::DIM,Tetra::NPI> _gradV;
calc_gradV(tet, _gradV);
gradV.push_back(_gradV);
Eigen::Matrix<double,Nodes::DIM,Tetra::NPI> _Hm;
calc_Hm(tet, _gradV, _Hm);
Hm.push_back(_Hm);
});
prepareExtras();
}
else
{
std::cerr << "Solver (STT) has not converged" << std::endl;
exit(1);
}
}
/** computes the gradient(V) for tetra tet */
void calc_gradV(Tetra::Tet const &tet, Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,Tetra::NPI>> _gradV) const
{
for (int npi = 0; npi < Tetra::NPI; npi++)
{
Eigen::Vector3d v(0,0,0);
for (int i = 0; i < Tetra::N; i++)
{
v += V[tet.ind[i]] * tet.da.row(i);
}
_gradV.col(npi) = v;
}
}
/** computes Hm contributions for each npi for tetrahedron tet */
void calc_Hm(Tetra::Tet const &tet, Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,Tetra::NPI>> _gradV,
Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,Tetra::NPI>> _Hm) const
{
Eigen::Matrix<double,Nodes::DIM,Tetra::NPI> p_g;
tet.getPtGauss(p_g);
for (int npi = 0; npi < Tetra::NPI; npi++)
{ _Hm.col(npi) = -p_stt.sigma * _gradV.col(npi).cross(p_g.col(npi)); }
}
private:
/** affect extraField function and extraCoeffs_BE function for all the tetrahedrons */
void prepareExtras(void)
{
std::for_each(
msh.tet.begin(), msh.tet.end(),
[this](Tetra::Tet &tet)
{
const int _idx = tet.idx;
tet.extraField = [this, _idx](Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,Tetra::NPI>> H)
{
for(int npi = 0; npi<Tetra::NPI; npi++) { H.col(npi) += Hm[_idx].col(npi); }
};
tet.extraCoeffs_BE = [this, &tet](double Js, Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,Tetra::NPI>> U,
Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,Tetra::NPI>> dUdx,
Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,Tetra::NPI>> dUdy,
Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,Tetra::NPI>> dUdz,
Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,Tetra::N>> BE)
{
const double prefactor = D0 / Nodes::sq(p_stt.lJ) / (gamma0*Js/mu0);
for (int npi = 0; npi < Tetra::NPI; npi++)
{
Eigen::Vector3d const &_gV = gradV[tet.idx].col(npi);
Eigen::Vector3d j_grad_u =
-p_stt.sigma
* Eigen::Vector3d(_gV.dot( Eigen::Vector3d(dUdx(Nodes::IDX_X,npi), dUdy(Nodes::IDX_X,npi), dUdz(Nodes::IDX_X,npi))),
_gV.dot( Eigen::Vector3d(dUdx(Nodes::IDX_Y,npi), dUdy(Nodes::IDX_Y,npi), dUdz(Nodes::IDX_Y,npi))),
_gV.dot( Eigen::Vector3d(dUdx(Nodes::IDX_Z,npi), dUdy(Nodes::IDX_Z,npi), dUdz(Nodes::IDX_Z,npi))));
Eigen::Vector3d m = pf * (ksi * j_grad_u + U.col(npi).cross(j_grad_u));
for (int i = 0; i < Tetra::N; i++)
{
const double ai_w = tet.weight[npi] * Tetra::a[i][npi];
BE.col(i) += ai_w*(Hm[tet.idx].col(npi) + prefactor*m);
}
} // end loop on npi
}; //end lambda
});//end for_each
}
/** ksi is in Thiaville notations beta_DW */
double ksi;
/** density of states */
double D0;
/** a prefactor for BE coefficient coefficients*/
double pf;
/** mesh object to store nodes, fac, tet, and others geometrical values related to the mesh (
* const ref ) */
Mesh::mesh msh;
/** spin transfer torque parameters */
STT p_stt;
/** if verbose set to true, some printing are sent to terminal */
const bool verbose;
/** maximum number of iteration for biconjugate stabilized gradient */
const unsigned int MAXITER; // fixed to 5000 in ref code
/** number of digits in the optional output file */
const int precision;
/** output .sol file name for electrostatic problem */
const std::string fileName;
/** electrostatic potential values for boundary conditions, V.size() is the size of the vector
* of nodes */
std::vector<double> V;
/** table of the gradients of the potential, gradV.size() is the number of tetra */
std::vector< Eigen::Matrix<double,Nodes::DIM,Tetra::NPI> > gradV;
/** table of the Hm vectors (contribution of the STT to the tet::integrales) ; Hm.size() is the
* number of tetra */
std::vector< Eigen::Matrix<double,Nodes::DIM,Tetra::NPI> > Hm;
/** basic informations on boundary conditions */
inline void infos(void)
{
std::cout << "sigma: " << p_stt.sigma << std::endl;
std::for_each(p_stt.boundaryCond.begin(), p_stt.boundaryCond.end(),
[](std::pair<std::string, double> const &p)
{ std::cout << "regName: " << p.first << "\tV :" << p.second << std::endl; });
}
/** fill matrix and vector to solve potential values on each node */
void prepareData(std::vector<Eigen::Triplet<double>> &Kw, Eigen::Ref<Eigen::VectorXd> Lw)
{
const double sigma = p_stt.sigma;
std::for_each(msh.tet.begin(), msh.tet.end(),
[this, sigma, &Kw](Tetra::Tet const &tet)
{
double K[Tetra::N][Tetra::N] = {{0}};
integrales(tet, sigma, K);
assembling_mat(tet, K, Kw);
});
double pot_val = 0; // we initialize pot_val to the average of the potentials set by the
// boundary conditions (does not seem to change convergence speed
// whatever value it is ..)
std::for_each(p_stt.boundaryCond.begin(), p_stt.boundaryCond.end(),
[&pot_val](auto const &it) { pot_val += it.second; });
pot_val /= p_stt.boundaryCond.size();
std::for_each(msh.fac.begin(), msh.fac.end(),
[this, pot_val, &Lw](Facette::Fac const &fac)
{
std::vector<double> L(Facette::N);
integrales(fac, pot_val, L);
assembling_vect(fac, L, Lw);
});
}
/** solver, using biconjugate stabilized gradient, with diagonal preconditionner and Dirichlet
* boundary conditions */
int solve(const double _tol)
{
const int NOD = msh.getNbNodes();
std::vector<Eigen::Triplet<double>> Kw;
Eigen::VectorXd Lw = Eigen::VectorXd::Zero(NOD);
prepareData(Kw, Lw);
Eigen::VectorXd Xw(NOD);
/* do something here with boundary conditions */
if (verbose)
{ std::cout << "line weighting..." << std::endl; }
Eigen::SparseMatrix<double> Kr(NOD,NOD);
Kr.setFromTriplets(Kw.begin(),Kw.end());
std::vector<double> maxRow(NOD,0);
//here we fill the vector maxRow with infinity norm : maxRow[i] = max{|Kr.row(i)|}
for(int i=0;i<NOD;i++)
{
for(int k=0;k<Kr.outerSize();++k)
for(Eigen::SparseMatrix<double>::InnerIterator it(Kr,k); it; ++it)
{ if((it.row() == i)&&(fabs(it.value()) > maxRow[i])) { maxRow[i] = fabs(it.value()); } }
}
for (int i = 0; i < NOD; i++)
{
double norme = maxRow[i];
Lw[i] /= norme;
Kr.row(i) /= norme;
}
if (verbose)
{ std::cout << "solving ..." << std::endl; }
Eigen::BiCGSTAB<Eigen::SparseMatrix<double>> _solver;
_solver.setTolerance(_tol);
_solver.setMaxIterations(MAXITER);
_solver.compute(Kr);
Eigen::VectorXd sol = _solver.solve(Lw);
for (int i=0;i<NOD;i++)
{ V[i]= sol(i); }
return (_solver.iterations() < MAXITER);
}
}; // end class electrostatSolver