Skip to content

Commit c35b437

Browse files
author
Roberto Di Remigio
committed
Form T^-1R explicitly
1 parent 4b142ea commit c35b437

File tree

3 files changed

+153
-28
lines changed

3 files changed

+153
-28
lines changed

‎src/solver/IEFSolver.cpp‎

Lines changed: 7 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -50,34 +50,28 @@ void IEFSolver::buildSystemMatrix_impl(const Cavity & cavity, const IGreensFunct
5050

5151
void IEFSolver::buildAnisotropicMatrix(const Cavity & cav, const IGreensFunction & gf_i, const IGreensFunction & gf_o)
5252
{
53-
Tepsilon_ = anisotropicTEpsilon(cav, gf_i, gf_o);
53+
TinvR_ = anisotropicIEFMatrix(cav, gf_i, gf_o);
5454
// Pack into a block diagonal matrix
5555
// The number of irreps in the group
5656
int nrBlocks = cav.pointGroup().nrIrrep();
5757
// The size of the irreducible portion of the cavity
5858
int dimBlock = cav.irreducible_size();
5959
// For the moment just packs into a std::vector<Eigen::MatrixXd>
60-
symmetryPacking(blockTepsilon_, Tepsilon_, dimBlock, nrBlocks);
61-
62-
Rinfinity_ = anisotropicRinfinity(cav, gf_i, gf_o);
63-
symmetryPacking(blockRinfinity_, Rinfinity_, dimBlock, nrBlocks);
60+
symmetryPacking(blockTinvR_, TinvR_, dimBlock, nrBlocks);
6461

6562
built_ = true;
6663
}
6764

6865
void IEFSolver::buildIsotropicMatrix(const Cavity & cav, const IGreensFunction & gf_i, const IGreensFunction & gf_o)
6966
{
70-
Tepsilon_ = isotropicTEpsilon(cav, gf_i, profiles::epsilon(gf_o.permittivity()));
67+
TinvR_ = isotropicIEFMatrix(cav, gf_i, profiles::epsilon(gf_o.permittivity()));
7168
// Pack into a block diagonal matrix
7269
// The number of irreps in the group
7370
int nrBlocks = cav.pointGroup().nrIrrep();
7471
// The size of the irreducible portion of the cavity
7572
int dimBlock = cav.irreducible_size();
7673
// For the moment just packs into a std::vector<Eigen::MatrixXd>
77-
symmetryPacking(blockTepsilon_, Tepsilon_, dimBlock, nrBlocks);
78-
79-
Rinfinity_ = isotropicRinfinity(cav, gf_i);
80-
symmetryPacking(blockRinfinity_, Rinfinity_, dimBlock, nrBlocks);
74+
symmetryPacking(blockTinvR_, TinvR_, dimBlock, nrBlocks);
8175

8276
built_ = true;
8377
}
@@ -87,12 +81,12 @@ Eigen::VectorXd IEFSolver::computeCharge_impl(const Eigen::VectorXd & potential,
8781
// The potential and charge vector are of dimension equal to the
8882
// full dimension of the cavity. We have to select just the part
8983
// relative to the irrep needed.
90-
int fullDim = Rinfinity_.rows();
84+
int fullDim = TinvR_.rows();
9185
Eigen::VectorXd charge = Eigen::VectorXd::Zero(fullDim);
92-
int nrBlocks = blockRinfinity_.size();
86+
int nrBlocks = blockTinvR_.size();
9387
int irrDim = fullDim/nrBlocks;
9488
charge.segment(irrep*irrDim, irrDim) =
95-
- blockTepsilon_[irrep].partialPivLu().solve(blockRinfinity_[irrep] * potential.segment(irrep*irrDim, irrDim));
89+
- blockTinvR_[irrep] * potential.segment(irrep*irrDim, irrDim);
9690
// Symmetrize ASC charge := (charge + charge*)/2
9791
if (hermitivitize_) hermitivitize(charge);
9892

‎src/solver/IEFSolver.hpp‎

Lines changed: 5 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -45,11 +45,7 @@ class IGreensFunction;
4545
* \author Luca Frediani, Roberto Di Remigio
4646
* \date 2011, 2015, 2016
4747
*
48-
* \note We store the unsymmetrized T(epsilon) and Rinfinity matrices.
49-
* The ASC is obtained by multiplying the MEP by Rinfinity and then using a partially
50-
* pivoted LU decomposition of T(epsilon) on the resulting vector.
51-
* The ASC is then symmetrized. This avoids computing and storing the inverse
52-
* explicitly.
48+
* \note We store the unsymmetrized T^-1R matrix.
5349
*/
5450

5551
class IEFSolver : public PCMSolver
@@ -79,14 +75,10 @@ class IEFSolver : public PCMSolver
7975
private:
8076
/*! Whether the system matrix has to be symmetrized */
8177
bool hermitivitize_;
82-
/*! T(epsilon) matrix, not symmetry blocked */
83-
Eigen::MatrixXd Tepsilon_;
84-
/*! T(epsilon) matrix, symmetry blocked form */
85-
std::vector<Eigen::MatrixXd> blockTepsilon_;
86-
/*! R_infinity matrix, not symmetry blocked */
87-
Eigen::MatrixXd Rinfinity_;
88-
/*! R_infinity matrix, symmetry blocked form */
89-
std::vector<Eigen::MatrixXd> blockRinfinity_;
78+
/*! T^-1R matrix, not symmetry blocked */
79+
Eigen::MatrixXd TinvR_;
80+
/*! T^-1R matrix, symmetry blocked form */
81+
std::vector<Eigen::MatrixXd> blockTinvR_;
9082

9183
/*! \brief Calculation of the PCM matrix
9284
* \param[in] cavity the cavity to be used

‎src/solver/SolverImpl.hpp‎

Lines changed: 141 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,146 @@
4141
* \date 2015
4242
*/
4343

44+
/*! \brief Builds the **anisotropic** IEFPCM matrix
45+
* \param[in] cav the discretized cavity
46+
* \param[in] gf_i Green's function inside the cavity
47+
* \param[in] gf_o Green's function outside the cavity
48+
* \return the \f$ \mathbf{K} = \mathbf{T}^{-1}\mathbf{R}\mathbf{A} \f$ matrix
49+
*
50+
* This function calculates the PCM matrix. We use the following definitions:
51+
* \f[
52+
* \begin{align}
53+
* \mathbf{T} &=
54+
* \left(2\pi\mathbf{I} - \mathbf{D}_\mathrm{e}\mathbf{A}\right)\mathbf{S}_\mathrm{i}
55+
* +\mathbf{S}_\mathrm{e}\left(2\pi\mathbf{I} +
56+
* \mathbf{A}\mathbf{D}_\mathrm{i}^\dagger\right) \\
57+
* \mathbf{R} &=
58+
* \left(2\pi\mathbf{A}^{-1} - \mathbf{D}_\mathrm{e}\right) -
59+
* \mathbf{S}_\mathrm{e}\mathbf{S}^{-1}_\mathrm{i}\left(2\pi\mathbf{A}^{-1}-\mathbf{D}_\mathrm{i}\right)
60+
* \end{align}
61+
* \f]
62+
* The matrix is not symmetrized and is not symmetry packed.
63+
*/
64+
inline Eigen::MatrixXd anisotropicIEFMatrix(const Cavity & cav, const IGreensFunction & gf_i, const IGreensFunction & gf_o)
65+
{
66+
// The total size of the cavity
67+
PCMSolverIndex cavitySize = cav.size();
68+
// The number of irreps in the group
69+
int nrBlocks = cav.pointGroup().nrIrrep();
70+
// The size of the irreducible portion of the cavity
71+
int dimBlock = cav.irreducible_size();
72+
73+
// Compute SI, DI and SE, DE on the whole cavity, regardless of symmetry
74+
TIMER_ON("Computing SI");
75+
Eigen::MatrixXd SI = gf_i.singleLayer(cav.elements());
76+
TIMER_OFF("Computing SI");
77+
TIMER_ON("Computing DI");
78+
Eigen::MatrixXd DI = gf_i.doubleLayer(cav.elements());
79+
TIMER_OFF("Computing DI");
80+
TIMER_ON("Computing SE");
81+
Eigen::MatrixXd SE = gf_o.singleLayer(cav.elements());
82+
TIMER_OFF("Computing SE");
83+
TIMER_ON("Computing DE");
84+
Eigen::MatrixXd DE = gf_o.doubleLayer(cav.elements());
85+
TIMER_OFF("Computing DE");
86+
87+
// Perform symmetry blocking
88+
// If the group is C1 avoid symmetry blocking, we will just pack the fullPCMMatrix
89+
// into "block diagonal" when all other manipulations are done.
90+
if (cav.pointGroup().nrGenerators() != 0) {
91+
TIMER_ON("Symmetry blocking");
92+
symmetryBlocking(DI, cavitySize, dimBlock, nrBlocks);
93+
symmetryBlocking(SI, cavitySize, dimBlock, nrBlocks);
94+
symmetryBlocking(DE, cavitySize, dimBlock, nrBlocks);
95+
symmetryBlocking(SE, cavitySize, dimBlock, nrBlocks);
96+
TIMER_OFF("Symmetry blocking");
97+
}
98+
99+
Eigen::MatrixXd a = cav.elementArea().asDiagonal();
100+
Eigen::MatrixXd Id = Eigen::MatrixXd::Identity(cavitySize, cavitySize);
101+
102+
TIMER_ON("Assemble T matrix");
103+
Eigen::MatrixXd T = ((2 * M_PI * Id - DE * a) * SI + SE * (2 * M_PI * Id + a * DI.adjoint().eval()));
104+
TIMER_OFF("Assemble T matrix");
105+
106+
TIMER_ON("Assemble R matrix");
107+
Eigen::MatrixXd R = ((2 * M_PI * Id - DE * a) - SE * SI.ldlt().solve((2 * M_PI * Id - DI * a)));
108+
TIMER_OFF("Assemble R matrix");
109+
110+
TIMER_ON("Assemble T^-1R matrix");
111+
Eigen::MatrixXd fullPCMMatrix = T.partialPivLu().solve(R);
112+
TIMER_OFF("Assemble T^-1R matrix");
113+
114+
return fullPCMMatrix;
115+
}
116+
117+
/*! \brief Builds the **isotropic** IEFPCM matrix
118+
* \param[in] cav the discretized cavity
119+
* \param[in] gf_i Green's function inside the cavity
120+
* \param[in] epsilon permittivity outside the cavity
121+
* \return the \f$ \mathbf{K} = \mathbf{T}^{-1}\mathbf{R}\mathbf{A} \f$ matrix
122+
*
123+
* This function calculates the PCM matrix. We use the following definitions:
124+
* \f[
125+
* \begin{align}
126+
* \mathbf{T} &=
127+
* \left(2\pi\frac{\varepsilon+1}{\varepsilon-1}\mathbf{I} - \mathbf{D}_\mathrm{i}\mathbf{A}\right)\mathbf{S}_\mathrm{i} \\
128+
* \mathbf{R} &=
129+
* \left(2\pi\mathbf{A}^{-1} - \mathbf{D}_\mathrm{i}\right)
130+
* \end{align}
131+
* \f]
132+
* The matrix is not symmetrized and is not symmetry packed.
133+
*/
134+
inline Eigen::MatrixXd isotropicIEFMatrix(const Cavity & cav, const IGreensFunction & gf_i, double epsilon)
135+
{
136+
// The total size of the cavity
137+
PCMSolverIndex cavitySize = cav.size();
138+
// The number of irreps in the group
139+
int nrBlocks = cav.pointGroup().nrIrrep();
140+
// The size of the irreducible portion of the cavity
141+
int dimBlock = cav.irreducible_size();
142+
143+
// Compute SI and DI on the whole cavity, regardless of symmetry
144+
TIMER_ON("Computing SI");
145+
Eigen::MatrixXd SI = gf_i.singleLayer(cav.elements());
146+
TIMER_OFF("Computing SI");
147+
TIMER_ON("Computing DI");
148+
Eigen::MatrixXd DI = gf_i.doubleLayer(cav.elements());
149+
TIMER_OFF("Computing DI");
150+
151+
// Perform symmetry blocking
152+
// If the group is C1 avoid symmetry blocking, we will just pack the fullPCMMatrix
153+
// into "block diagonal" when all other manipulations are done.
154+
if (cav.pointGroup().nrGenerators() != 0) {
155+
TIMER_ON("Symmetry blocking");
156+
symmetryBlocking(DI, cavitySize, dimBlock, nrBlocks);
157+
symmetryBlocking(SI, cavitySize, dimBlock, nrBlocks);
158+
TIMER_OFF("Symmetry blocking");
159+
}
160+
161+
Eigen::MatrixXd a = cav.elementArea().asDiagonal();
162+
Eigen::MatrixXd Id = Eigen::MatrixXd::Identity(cavitySize, cavitySize);
163+
164+
// Tq = -Rv -> q = -(T^-1 * R)v = -Kv
165+
// T = (2 * M_PI * fact * aInv - DI) * a * SI; R = (2 * M_PI * aInv - DI)
166+
// fullPCMMatrix_ = K = T^-1 * R * a
167+
// 1. Form T
168+
double fact = (epsilon + 1.0)/(epsilon - 1.0);
169+
TIMER_ON("Assemble T matrix");
170+
Eigen::MatrixXd T = (2 * M_PI * fact * Id - DI * a) * SI;
171+
TIMER_OFF("Assemble T matrix");
172+
173+
TIMER_ON("Assemble R matrix");
174+
Eigen::MatrixXd R = (2 * M_PI * Id - DI * a);
175+
TIMER_OFF("Assemble R matrix");
176+
177+
TIMER_ON("Assemble T^-1R matrix");
178+
Eigen::MatrixXd fullPCMMatrix = T.partialPivLu().solve(R);
179+
TIMER_OFF("Assemble T^-1R matrix");
180+
181+
return fullPCMMatrix;
182+
}
183+
44184
/*! \brief Builds the **anisotropic** \f$ \mathbf{T}_\varepsilon \f$ matrix
45185
* \param[in] cav the discretized cavity
46186
* \param[in] gf_i Green's function inside the cavity
@@ -172,8 +312,7 @@ inline Eigen::MatrixXd anisotropicRinfinity(const Cavity & cav, const IGreensFun
172312
Eigen::MatrixXd Id = Eigen::MatrixXd::Identity(cavitySize, cavitySize);
173313

174314
// Form R
175-
Eigen::PartialPivLU<Eigen::MatrixXd> SI_LU(SI);
176-
return ((2 * M_PI * Id - DE * a) - SE * SI_LU.inverse() * (2 * M_PI * Id - DI * a));
315+
return ((2 * M_PI * Id - DE * a) - SE * SI.ldlt().solve((2 * M_PI * Id - DI * a)));
177316
}
178317

179318
/*! \brief Builds the **isotropic** \f$ \mathbf{R}_\infty \f$ matrix

0 commit comments

Comments
 (0)