Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <algorithm>
7 : #include <cstddef>
8 : #include <fstream>
9 : #include <string>
10 : #include <tuple>
11 : #include <vector>
12 :
13 : #include "DataStructures/DataBox/DataBox.hpp"
14 : #include "DataStructures/DynamicMatrix.hpp"
15 : #include "DataStructures/DynamicVector.hpp"
16 : #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
17 : #include "NumericalAlgorithms/LinearSolver/BuildMatrix.hpp"
18 : #include "NumericalAlgorithms/LinearSolver/LinearSolver.hpp"
19 : #include "Options/Auto.hpp"
20 : #include "Options/String.hpp"
21 : #include "Parallel/Tags/ArrayIndex.hpp"
22 : #include "Utilities/ErrorHandling/Error.hpp"
23 : #include "Utilities/Gsl.hpp"
24 : #include "Utilities/MakeWithValue.hpp"
25 : #include "Utilities/NoSuchType.hpp"
26 : #include "Utilities/Serialization/CharmPupable.hpp"
27 : #include "Utilities/Serialization/PupStlCpp17.hpp"
28 : #include "Utilities/TMPL.hpp"
29 :
30 : /// \cond
31 : template <size_t VolumeDim>
32 : class ElementId;
33 : /// \endcond
34 :
35 : namespace LinearSolver::Serial {
36 :
37 : /// \cond
38 : template <typename ValueType, typename LinearSolverRegistrars>
39 : struct ExplicitInverse;
40 : /// \endcond
41 :
42 1 : namespace Registrars {
43 : /// Registers the `LinearSolver::Serial::ExplicitInverse` linear solver
44 : template <typename ValueType>
45 1 : struct ExplicitInverse {
46 : template <typename LinearSolverRegistrars>
47 0 : using f = Serial::ExplicitInverse<ValueType, LinearSolverRegistrars>;
48 : };
49 : } // namespace Registrars
50 :
51 : /*!
52 : * \brief Linear solver that builds a matrix representation of the linear
53 : * operator and inverts it directly
54 : *
55 : * This solver first constructs an explicit matrix representation by "sniffing
56 : * out" the operator, i.e. feeding it with unit vectors, and then directly
57 : * inverts the matrix. The result is an operator that solves the linear problem
58 : * in a single step. This means that each element has a large initialization
59 : * cost, but all successive solves converge immediately.
60 : *
61 : * \par Advice on using this linear solver:
62 : * - This solver is entirely agnostic to the structure of the linear operator.
63 : * It is usually better to implement a linear solver that is specialized for
64 : * your linear operator to take advantage of its properties. For example, if
65 : * the operator has a tensor-product structure, the linear solver might take
66 : * advantage of that. Only use this solver if no alternatives are available
67 : * and if you have verified that it speeds up your solves.
68 : * - Since this linear solver stores the full inverse operator matrix it can
69 : * have significant memory demands. For example, an operator representing a 3D
70 : * first-order Elasticity system (9 variables) discretized on 12 grid points
71 : * per dimension requires ca. 2GB of memory (per element) to store the matrix,
72 : * scaling quadratically with the number of variables and with a power of 6
73 : * with the number of grid points per dimension. Therefore, make sure to
74 : * distribute the elements on a sufficient number of nodes to meet the memory
75 : * requirements.
76 : * - This linear solver can be `reset()` when the operator changes (e.g. in each
77 : * nonlinear-solver iteration). However, when using this solver as
78 : * preconditioner it can be advantageous to avoid the reset and the
79 : * corresponding cost of re-building the matrix and its inverse if the
80 : * operator only changes "a little". In that case the preconditioner solves
81 : * subdomain problems only approximately, but possibly still sufficiently to
82 : * provide effective preconditioning.
83 : */
84 : template <typename ValueType,
85 : typename LinearSolverRegistrars =
86 : tmpl::list<Registrars::ExplicitInverse<ValueType>>>
87 1 : class ExplicitInverse : public LinearSolver<LinearSolverRegistrars> {
88 : private:
89 0 : using Base = LinearSolver<LinearSolverRegistrars>;
90 :
91 : public:
92 0 : struct WriteMatrixToFile {
93 0 : using type = Options::Auto<std::string, Options::AutoLabel::None>;
94 0 : static constexpr Options::String help =
95 : "Write the matrix representation of the linear operator to a "
96 : "space-delimited CSV file with this name. A '.txt' extension will be "
97 : "added. Also a suffix with the element ID will be added if this linear "
98 : "solver runs on an array element, so one file per element will be "
99 : "written.";
100 : };
101 0 : struct ComplexShift {
102 0 : using type = double;
103 0 : static constexpr Options::String help =
104 : "Subtracts i*ComplexShift from the operator. "
105 : "This can help dampen high-frequency modes when preconditioning "
106 : "complex-valued problems.";
107 : };
108 :
109 0 : using options = tmpl::flatten<tmpl::list<
110 : WriteMatrixToFile,
111 : tmpl::conditional_t<std::is_same_v<ValueType, std::complex<double>>,
112 : ComplexShift, tmpl::list<>>>>;
113 0 : static constexpr Options::String help =
114 : "Build a matrix representation of the linear operator and invert it "
115 : "directly. This means that the first solve has a large initialization "
116 : "cost, but all subsequent solves converge immediately.";
117 :
118 0 : ExplicitInverse(const ExplicitInverse& /*rhs*/) = default;
119 0 : ExplicitInverse& operator=(const ExplicitInverse& /*rhs*/) = default;
120 0 : ExplicitInverse(ExplicitInverse&& /*rhs*/) = default;
121 0 : ExplicitInverse& operator=(ExplicitInverse&& /*rhs*/) = default;
122 0 : ~ExplicitInverse() = default;
123 :
124 0 : explicit ExplicitInverse(
125 : std::optional<std::string> matrix_filename = std::nullopt,
126 : const double complex_shift = 0.0)
127 : : matrix_filename_(std::move(matrix_filename)),
128 : complex_shift_(complex_shift) {}
129 :
130 : /// \cond
131 : explicit ExplicitInverse(CkMigrateMessage* m) : Base(m) {}
132 : using PUP::able::register_constructor;
133 : WRAPPED_PUPable_decl_template(ExplicitInverse); // NOLINT
134 : /// \endcond
135 :
136 : /*!
137 : * \brief Solve the equation \f$Ax=b\f$ by explicitly constructing the
138 : * operator matrix \f$A\f$ and its inverse. The first solve is computationally
139 : * expensive and successive solves are cheap.
140 : *
141 : * Building a matrix representation of the `linear_operator` requires
142 : * iterating over the `SourceType` (which is also the type returned by the
143 : * `linear_operator`) in a consistent way. This can be non-trivial for
144 : * heterogeneous data structures because it requires they define a data
145 : * ordering. Specifically, the `SourceType` must have a `size()` function as
146 : * well as `begin()` and `end()` iterators that point into the data. If the
147 : * iterators have a `reset()` function it is used to avoid repeatedly
148 : * re-creating the `begin()` iterator. The `reset()` function must not
149 : * invalidate the `end()` iterator.
150 : */
151 : template <typename LinearOperator, typename VarsType, typename SourceType,
152 : typename... OperatorArgs>
153 1 : Convergence::HasConverged solve(
154 : gsl::not_null<VarsType*> solution, const LinearOperator& linear_operator,
155 : const SourceType& source,
156 : const std::tuple<OperatorArgs...>& operator_args = std::tuple{}) const;
157 :
158 : /// Flags the operator to require re-initialization. No memory is released.
159 : /// Call this function to rebuild the solver when the operator changed.
160 1 : void reset() override { size_ = std::numeric_limits<size_t>::max(); }
161 :
162 : /// Size of the operator. The stored matrix will have `size^2` entries.
163 1 : size_t size() const { return size_; }
164 :
165 : /// The matrix representation of the solver. This matrix approximates the
166 : /// inverse of the subdomain operator.
167 : const blaze::DynamicMatrix<ValueType, blaze::columnMajor>&
168 1 : matrix_representation() const {
169 : return inverse_;
170 : }
171 :
172 : // NOLINTNEXTLINE(google-runtime-references)
173 0 : void pup(PUP::er& p) override {
174 : p | matrix_filename_;
175 : p | complex_shift_;
176 : p | size_;
177 : p | inverse_;
178 : if (p.isUnpacking() and size_ != std::numeric_limits<size_t>::max()) {
179 : source_workspace_.resize(size_);
180 : solution_workspace_.resize(size_);
181 : }
182 : }
183 :
184 0 : std::unique_ptr<Base> get_clone() const override {
185 : return std::make_unique<ExplicitInverse>(*this);
186 : }
187 :
188 : private:
189 0 : std::optional<std::string> matrix_filename_{};
190 0 : double complex_shift_ = std::numeric_limits<double>::signaling_NaN();
191 : // Caches for successive solves of the same operator
192 : // NOLINTNEXTLINE(spectre-mutable)
193 0 : mutable size_t size_ = std::numeric_limits<size_t>::max();
194 : // We currently store the matrix representation in a dense matrix because
195 : // Blaze doesn't support the inversion of sparse matrices (yet).
196 : // NOLINTNEXTLINE(spectre-mutable)
197 0 : mutable blaze::DynamicMatrix<ValueType, blaze::columnMajor> inverse_{};
198 :
199 : // Buffers to avoid re-allocating memory for applying the operator
200 : // NOLINTNEXTLINE(spectre-mutable)
201 0 : mutable blaze::DynamicVector<ValueType> source_workspace_{};
202 : // NOLINTNEXTLINE(spectre-mutable)
203 0 : mutable blaze::DynamicVector<ValueType> solution_workspace_{};
204 : };
205 :
206 : template <typename ValueType, typename LinearSolverRegistrars>
207 : template <typename LinearOperator, typename VarsType, typename SourceType,
208 : typename... OperatorArgs>
209 : Convergence::HasConverged
210 0 : ExplicitInverse<ValueType, LinearSolverRegistrars>::solve(
211 : const gsl::not_null<VarsType*> solution,
212 : const LinearOperator& linear_operator, const SourceType& source,
213 : const std::tuple<OperatorArgs...>& operator_args) const {
214 : if (UNLIKELY(size_ == std::numeric_limits<size_t>::max())) {
215 : const auto& used_for_size = source;
216 : size_ = used_for_size.size();
217 : source_workspace_.resize(size_);
218 : solution_workspace_.resize(size_);
219 : inverse_.resize(size_, size_);
220 : // Construct explicit matrix representation by "sniffing out" the operator,
221 : // i.e. feeding it unit vectors
222 : auto operand_buffer = make_with_value<VarsType>(used_for_size, 0.);
223 : auto result_buffer = make_with_value<SourceType>(used_for_size, 0.);
224 : build_matrix(make_not_null(&inverse_), make_not_null(&operand_buffer),
225 : make_not_null(&result_buffer), linear_operator, operator_args);
226 : // Write to file before inverting
227 : if (UNLIKELY(matrix_filename_.has_value())) {
228 : std::string filename_suffix{};
229 : using DataBoxType =
230 : std::decay_t<tmpl::front<tmpl::list<OperatorArgs..., NoSuchType>>>;
231 : if constexpr (tt::is_a_v<db::DataBox, DataBoxType>) {
232 : tmpl::for_each<tmpl::range<size_t, 1, 4>>(
233 : [&]<size_t Dim>(tmpl::type_<tmpl::size_t<Dim>> /*meta*/) {
234 : using index_tag = Parallel::Tags::ArrayIndex<ElementId<Dim>>;
235 : if constexpr (db::tag_is_retrievable_v<index_tag, DataBoxType>) {
236 : const auto& box = std::get<0>(operator_args);
237 : filename_suffix = "_" + get_output(db::get<index_tag>(box));
238 : } else {
239 : (void)operator_args;
240 : }
241 : });
242 : } else {
243 : (void)operator_args;
244 : }
245 : std::ofstream matrix_file(matrix_filename_.value() + filename_suffix +
246 : ".txt");
247 : write_csv(matrix_file, inverse_, " ");
248 : }
249 : // Apply complex shift
250 : if constexpr (std::is_same_v<ValueType, std::complex<double>>) {
251 : for (size_t i = 0; i < size_; ++i) {
252 : inverse_(i, i) -= std::complex<double>(0.0, complex_shift_);
253 : }
254 : }
255 : // Directly invert the matrix
256 : try {
257 : blaze::invert(inverse_);
258 : } catch (const std::exception& e) {
259 : ERROR("Could not invert subdomain matrix (size " << size_
260 : << "): " << e.what());
261 : }
262 : }
263 : // Copy source into contiguous workspace. In cases where the source and
264 : // solution data are already stored contiguously we might avoid the copy and
265 : // the associated workspace memory. However, compared to the cost of building
266 : // and storing the matrix this is likely insignificant.
267 : std::copy(source.begin(), source.end(), source_workspace_.begin());
268 : // Apply inverse
269 : solution_workspace_ = inverse_ * source_workspace_;
270 : // Reconstruct solution data from contiguous workspace
271 : std::copy(solution_workspace_.begin(), solution_workspace_.end(),
272 : solution->begin());
273 : return {0, 0};
274 : }
275 :
276 : /// \cond
277 : // NOLINTBEGIN
278 : template <typename ValueType, typename LinearSolverRegistrars>
279 : PUP::able::PUP_ID
280 : ExplicitInverse<ValueType, LinearSolverRegistrars>::my_PUP_ID = 0;
281 : // NOLINTEND
282 : /// \endcond
283 :
284 : } // namespace LinearSolver::Serial
|
