-
Notifications
You must be signed in to change notification settings - Fork 34
/
kerkerSolverProblemDevice.h
246 lines (194 loc) · 7.02 KB
/
kerkerSolverProblemDevice.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
// ---------------------------------------------------------------------
//
// Copyright (c) 2019-2020 The Regents of the University of Michigan and DFT-FE
// authors.
//
// This file is part of the DFT-FE code.
//
// The DFT-FE code is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the DFT-FE distribution.
//
// ---------------------------------------------------------------------
//
#if defined(DFTFE_WITH_DEVICE)
# ifndef kerkerSolverProblemDevice_H_
# define kerkerSolverProblemDevice_H_
# include <linearSolverProblemDevice.h>
# include <triangulationManager.h>
# include <constraintMatrixInfo.h>
# include <MemoryStorage.h>
# include <dftUtils.h>
# include <FEBasisOperations.h>
# include <MatrixFree.h>
namespace dftfe
{
/**
* @brief helmholtz solver problem class template. template parameter FEOrderElectro
* is the finite element polynomial order for electrostatics
*
* @author Gourab Panigrahi
*/
template <unsigned int FEOrderElectro>
class kerkerSolverProblemDevice : public linearSolverProblemDevice
{
public:
/// Constructor
kerkerSolverProblemDevice(const MPI_Comm &mpi_comm_parent,
const MPI_Comm &mpi_comm_domain);
/**
* @brief initialize the matrix-free data structures
*
* @param matrixFreeData structure to hold quadrature rule, constraints vector and appropriate dofHandler
* @param constraintMatrix to hold constraints in the given problem
* @param x vector to be initialized using matrix-free object
*
*/
void
init(
std::shared_ptr<
dftfe::basis::
FEBasisOperations<double, double, dftfe::utils::MemorySpace::DEVICE>>
& basisOperationsPtr,
dealii::AffineConstraints<double> &constraintMatrix,
distributedCPUVec<double> & x,
double kerkerMixingParameter,
const unsigned int matrixFreeVectorComponent,
const unsigned int matrixFreeQuadratureComponent);
/**
* @brief reinitialize data structures .
*
* @param x vector to store initial guess and solution
* @param gradResidualValues stores the gradient of difference of input electron-density and output electron-density
* @param kerkerMixingParameter used in Kerker mixing scheme which usually represents Thomas Fermi wavevector (k_{TF}**2).
*
*/
void
reinit(
distributedCPUVec<double> &x,
const dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>
&quadPointValues);
/**
* @brief get the reference to x field
*
* @return reference to x field. Assumes x field data structure is already initialized
*/
distributedDeviceVec<double> &
getX();
/**
* @brief get the reference to Preconditioner
*
* @return reference to Preconditioner
*/
distributedDeviceVec<double> &
getPreconditioner();
/**
* @brief Compute A matrix multipled by x.
*
*/
void
computeAX(distributedDeviceVec<double> &Ax,
distributedDeviceVec<double> &x);
void
setX();
/**
* @brief Compute right hand side vector for the problem Ax = rhs.
*
* @param rhs vector for the right hand side values
*/
void
computeRhs(distributedCPUVec<double> &rhs);
/**
* @brief distribute x to the constrained nodes.
*
*/
void
distributeX();
/**
* @brief Copies x from Device to Host
*
*/
void
copyXfromDeviceToHost();
/// function needed by dealii to mimic SparseMatrix for Jacobi
/// preconditioning
void
subscribe(std::atomic<bool> *const validity,
const std::string & identifier = "") const {};
/// function needed by dealii to mimic SparseMatrix for Jacobi
/// preconditioning
void
unsubscribe(std::atomic<bool> *const validity,
const std::string & identifier = "") const {};
/// function needed by dealii to mimic SparseMatrix
bool
operator!=(double val) const
{
return true;
};
private:
/**
* @brief Sets up the constraints matrix
*
*/
void
setupConstraints();
/**
* @brief Compute the diagonal of A.
*
*/
void
computeDiagonalA();
/// storage for diagonal of the A matrix
distributedCPUVec<double> d_diagonalA;
distributedDeviceVec<double> d_diagonalAdevice;
/// pointer to the x vector being solved for
distributedCPUVec<double> *d_xPtr;
/// Device x vector being solved for
distributedDeviceVec<double> d_xDevice;
// number of cells local to each mpi task, number of degrees of freedom
// locally owned and total degrees of freedom including ghost
int d_nLocalCells, d_xLocalDof, d_xLen;
// kerker mixing constant
double d_gamma;
/// pointer to MatrixFree Base object
std::unique_ptr<dftfe::MatrixFreeBase> d_matrixFreeBasePtr;
// shape function value, gradient, jacobian and map for matrixfree
dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::DEVICE>
d_shapeFunction, d_jacobianFactor;
dftfe::utils::MemoryStorage<int, dftfe::utils::MemorySpace::DEVICE> d_map;
// Pointers to shape function value, gradient, jacobian and map for
// matrixfree
double *d_shapeFunctionPtr;
double *d_jacobianFactorPtr;
int * d_mapPtr;
// constraints
dftUtils::constraintMatrixInfo<dftfe::utils::MemorySpace::DEVICE>
d_constraintsTotalPotentialInfo;
/// matrix free index required to access the DofHandler and
/// dealii::AffineConstraints<double> objects corresponding to the problem
unsigned int d_matrixFreeVectorComponent;
/// matrix free quadrature index
unsigned int d_matrixFreeQuadratureComponent;
/// pointer to electron density cell and grad residual data
const dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>
* d_residualQuadValuesPtr;
const dealii::DoFHandler<3> * d_dofHandlerPRefinedPtr;
const dealii::AffineConstraints<double> *d_constraintMatrixPRefinedPtr;
const dealii::MatrixFree<3, double> * d_matrixFreeDataPRefinedPtr;
std::shared_ptr<
dftfe::basis::
FEBasisOperations<double, double, dftfe::utils::MemorySpace::DEVICE>>
d_basisOperationsPtr;
const MPI_Comm d_mpiCommParent;
const MPI_Comm mpi_communicator;
const unsigned int n_mpi_processes;
const unsigned int this_mpi_process;
dealii::ConditionalOStream pcout;
};
} // namespace dftfe
# endif // kerkerSolverProblemDevice_H_
#endif