Line data Source code
1 : // SPDX-FileCopyrightText: 2025 Pairinteraction Developers
2 : // SPDX-License-Identifier: LGPL-3.0-or-later
3 :
4 : #include "pairinteraction/system/GreenTensor.hpp"
5 :
6 : #include "pairinteraction/utils/eigen_assertion.hpp"
7 : #include "pairinteraction/utils/eigen_compat.hpp"
8 : #include "pairinteraction/utils/spherical.hpp"
9 : #include "pairinteraction/utils/traits.hpp"
10 :
11 : #include <Eigen/Dense>
12 : #include <Eigen/Sparse>
13 : #include <complex>
14 : #include <map>
15 : #include <spdlog/spdlog.h>
16 : #include <unsupported/Eigen/Splines>
17 :
18 : namespace pairinteraction {
19 : template <typename Scalar>
20 24 : GreenTensor<Scalar>::ConstantEntry::ConstantEntry(int row, int col, Scalar val)
21 24 : : row_(row), col_(col), val_(val) {}
22 :
23 : template <typename Scalar>
24 27 : Scalar GreenTensor<Scalar>::ConstantEntry::val() const {
25 27 : return val_;
26 : }
27 :
28 : template <typename Scalar>
29 54 : int GreenTensor<Scalar>::ConstantEntry::row() const noexcept {
30 54 : return row_;
31 : }
32 :
33 : template <typename Scalar>
34 54 : int GreenTensor<Scalar>::ConstantEntry::col() const noexcept {
35 54 : return col_;
36 : }
37 :
38 : template <typename Scalar>
39 0 : GreenTensor<Scalar>::OmegaDependentEntry::OmegaDependentEntry(int row, int col,
40 : Eigen::Spline<real_t, 1> real_spline,
41 : Eigen::Spline<real_t, 1> imag_spline)
42 0 : : row_(row), col_(col), real_spline(std::move(real_spline)),
43 0 : imag_spline(std::move(imag_spline)) {}
44 :
45 : template <typename Scalar>
46 0 : Scalar GreenTensor<Scalar>::OmegaDependentEntry::val(double omega) const {
47 : if constexpr (traits::NumTraits<Scalar>::is_complex_v) {
48 0 : return {real_spline(omega)(0), imag_spline(omega)(0)};
49 : } else {
50 0 : return real_spline(omega)(0);
51 : }
52 : }
53 :
54 : template <typename Scalar>
55 0 : int GreenTensor<Scalar>::OmegaDependentEntry::row() const noexcept {
56 0 : return row_;
57 : }
58 :
59 : template <typename Scalar>
60 0 : int GreenTensor<Scalar>::OmegaDependentEntry::col() const noexcept {
61 0 : return col_;
62 : }
63 :
64 : template <typename Scalar>
65 8 : void GreenTensor<Scalar>::create_entries_from_cartesian(
66 : int kappa1, int kappa2, const Eigen::MatrixX<Scalar> &tensor_in_cartesian_coordinates) {
67 :
68 8 : const real_t scale = tensor_in_cartesian_coordinates.norm();
69 9 : const real_t numerical_precision = 100 * scale * std::numeric_limits<real_t>::epsilon();
70 :
71 18 : Eigen::SparseMatrix<complex_t> tensor =
72 8 : (spherical::get_transformator<complex_t>(kappa1) * tensor_in_cartesian_coordinates *
73 8 : spherical::get_transformator<complex_t>(kappa2).adjoint())
74 9 : .sparseView(1, numerical_precision);
75 :
76 8 : std::vector<Entry> entries;
77 34 : for (int k = 0; k < tensor.outerSize(); ++k) {
78 50 : for (typename Eigen::SparseMatrix<complex_t>::InnerIterator it(tensor, k); it; ++it) {
79 : if constexpr (traits::NumTraits<Scalar>::is_complex_v) {
80 0 : entries.emplace_back(ConstantEntry(it.row(), it.col(), it.value()));
81 : } else {
82 23 : entries.emplace_back(ConstantEntry(it.row(), it.col(), it.value().real()));
83 27 : assert(abs(it.value().imag()) < numerical_precision);
84 : }
85 : }
86 : }
87 9 : entries_map[{kappa1, kappa2}] = std::move(entries);
88 8 : }
89 :
90 : template <typename Scalar>
91 0 : void GreenTensor<Scalar>::create_entries_from_cartesian(
92 : int kappa1, int kappa2,
93 : const std::vector<Eigen::MatrixX<Scalar>> &tensors_in_cartesian_coordinates,
94 : const std::vector<double> &omegas) {
95 :
96 0 : if (tensors_in_cartesian_coordinates.size() != omegas.size()) {
97 0 : throw std::invalid_argument("The number of tensors and omegas must match.");
98 : }
99 :
100 0 : if (tensors_in_cartesian_coordinates.size() < 4) {
101 0 : throw std::invalid_argument(
102 : "At least 4 tensors are required for the applied cubic spline interpolation.");
103 : }
104 :
105 0 : auto num_knots = static_cast<int>(omegas.size());
106 0 : Eigen::Map<const Eigen::RowVectorXd> knots(omegas.data(), num_knots);
107 :
108 0 : constexpr int spline_degree = 3; // cubic spline interpolation
109 :
110 : // Temporary storage wih key = (row, col) and value = vector of one double per omega
111 0 : std::map<std::pair<int, int>, std::pair<Eigen::RowVectorXd, Eigen::RowVectorXd>> temp_map;
112 0 : for (int idx = 0; idx < num_knots; ++idx) {
113 :
114 0 : const real_t scale = tensors_in_cartesian_coordinates[idx].norm();
115 0 : const real_t numerical_precision = 100 * scale * std::numeric_limits<real_t>::epsilon();
116 :
117 0 : Eigen::SparseMatrix<complex_t> tensor =
118 0 : (spherical::get_transformator<complex_t>(kappa1) *
119 0 : tensors_in_cartesian_coordinates[idx] *
120 0 : spherical::get_transformator<complex_t>(kappa2).adjoint())
121 0 : .sparseView(1, numerical_precision);
122 :
123 0 : for (int k = 0; k < tensor.outerSize(); ++k) {
124 0 : for (typename Eigen::SparseMatrix<complex_t>::InnerIterator it(tensor, k); it; ++it) {
125 0 : std::pair<int, int> key{it.row(), it.col()};
126 0 : auto &[vec_real, vec_imag] =
127 : temp_map
128 0 : .try_emplace(key, Eigen::RowVectorXd::Zero(num_knots),
129 : Eigen::RowVectorXd::Zero(num_knots))
130 0 : .first->second;
131 0 : vec_real(idx) = it.value().real();
132 : if constexpr (traits::NumTraits<Scalar>::is_complex_v) {
133 0 : vec_imag(idx) = it.value().imag();
134 : } else {
135 0 : assert(abs(it.value().imag()) < numerical_precision);
136 : }
137 : }
138 : }
139 : }
140 :
141 : // Set the green tensor entries with spline interpolation
142 0 : std::vector<Entry> entries;
143 0 : entries.reserve(temp_map.size());
144 0 : for (const auto &[key, value] : temp_map) {
145 0 : const auto &[vec_real, vec_imag] = value;
146 0 : const auto &[row, col] = key;
147 :
148 0 : Eigen::Spline<real_t, 1> real_spline =
149 : Eigen::SplineFitting<Eigen::Spline<real_t, 1>>::Interpolate(vec_real, spline_degree,
150 : knots);
151 :
152 0 : Eigen::Spline<real_t, 1> imag_spline;
153 : if constexpr (traits::NumTraits<Scalar>::is_complex_v) {
154 0 : imag_spline = Eigen::SplineFitting<Eigen::Spline<real_t, 1>>::Interpolate(
155 : vec_imag, spline_degree, knots);
156 : }
157 :
158 0 : entries.emplace_back(
159 0 : OmegaDependentEntry(row, col, std::move(real_spline), std::move(imag_spline)));
160 : }
161 0 : entries_map[{kappa1, kappa2}] = std::move(entries);
162 0 : }
163 :
164 : template <typename Scalar>
165 : const std::vector<typename GreenTensor<Scalar>::Entry> &
166 88 : GreenTensor<Scalar>::get_spherical_entries(int kappa1, int kappa2) const {
167 88 : if (auto it = entries_map.find({kappa1, kappa2}); it != entries_map.end()) {
168 25 : return it->second;
169 : }
170 63 : static const std::vector<Entry> empty_entries;
171 63 : return empty_entries;
172 : }
173 :
174 : // Explicit instantiations
175 : template class GreenTensor<double>;
176 : template class GreenTensor<std::complex<double>>;
177 : } // namespace pairinteraction
|
