Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #pragma once
11 :
12 : #include "MooseTypes.h"
13 : #include "MooseUtils.h"
14 : #include "DataIO.h"
15 : #include "libmesh/int_range.h"
16 : #include "libmesh/nanoflann.hpp"
17 : #include <cstdlib>
18 : #include <fstream>
19 : #include <iostream>
20 : #include <memory>
21 : #include <vector>
22 : #include <tuple>
23 : #include <optional>
24 : #include <exception>
25 :
26 : template <typename T>
27 : class ValueCache;
28 : template <typename T>
29 : void dataStore(std::ostream & stream, ValueCache<T> & c, void * context);
30 : template <typename T>
31 : void dataLoad(std::istream & stream, ValueCache<T> & c, void * context);
32 :
33 : /**
34 : * ValueCache is a generic helper template to implement an unstructured data
35 : * cache, where arbitrary result types can be placed in an n-dimensional space of
36 : * real numbers. Upon lookup the closest result to a given point in n-space is
37 : * returned. Applications include caching reasonable initial guesses for local
38 : * Newton solves as found in certain phase field models as well as thermodynamic
39 : * Gibbs energy minimization.
40 : */
41 : template <typename T>
42 : class ValueCache
43 : {
44 : public:
45 : /// Construct a ValueCache with indices in in_dim dimensions
46 : ValueCache(std::size_t in_dim, std::size_t max_leaf_size = 10);
47 :
48 : /// Same as above, but provide a filename for storing/restoring the cache content between runs
49 : ValueCache(const std::string & file_name, std::size_t in_dim, std::size_t max_leaf_size = 10);
50 :
51 : /// Object destructor. If the cache was constructed with a file name, it gets written here.
52 : ~ValueCache();
53 :
54 : /// insert a new value out_value at the position in_val
55 : void insert(const std::vector<Real> & in_val, const T & out_val);
56 :
57 : /// get a single neighbor of in_val along with stored values and distances
58 : std::tuple<const std::vector<Real> &, const T &, Real>
59 : getNeighbor(const std::vector<Real> & in_val);
60 :
61 : /// get a list of up to k neighbors of in_val along with stored values and distances
62 : std::vector<std::tuple<const std::vector<Real> &, const T &, Real>>
63 : getNeighbors(const std::vector<Real> & in_val, const std::size_t k);
64 :
65 : /// return the number of cache entries
66 : std::size_t size();
67 :
68 : /// remove all data from the cache
69 : void clear();
70 :
71 : ///@{ Nanoflann interface functions
72 : std::size_t kdtree_get_point_count() const;
73 : Real kdtree_get_pt(const std::size_t idx, const std::size_t dim) const;
74 : template <class BBOX>
75 : bool kdtree_get_bbox(BBOX & bb) const;
76 : ///@}
77 :
78 : protected:
79 : /// rebuild the kd-tree from scratch and update the bounding box
80 : void rebuildTree();
81 :
82 : using KdTreeT =
83 : nanoflann::KDTreeSingleIndexDynamicAdaptor<nanoflann::L2_Simple_Adaptor<Real, ValueCache<T>>,
84 : ValueCache<T>,
85 : -1,
86 : std::size_t>;
87 :
88 : std::vector<std::pair<std::vector<Real>, T>> _location_data;
89 : std::unique_ptr<KdTreeT> _kd_tree;
90 :
91 : const std::size_t _in_dim;
92 : const std::size_t _max_leaf_size;
93 : const std::size_t _max_subtrees;
94 :
95 : /// file name for persistent store/restore of the cache
96 : std::optional<std::string> _persistent_storage_file;
97 :
98 : /// bounding box (updated upon insertion)
99 : std::vector<std::pair<Real, Real>> _bbox;
100 :
101 : friend void dataStore<T>(std::ostream & stream, ValueCache<T> & c, void * context);
102 : friend void dataLoad<T>(std::istream & stream, ValueCache<T> & c, void * context);
103 : };
104 :
105 : template <typename T>
106 6 : ValueCache<T>::ValueCache(std::size_t in_dim, std::size_t max_leaf_size)
107 6 : : _kd_tree(nullptr), _in_dim(in_dim), _max_leaf_size(max_leaf_size), _max_subtrees(100)
108 : {
109 6 : }
110 :
111 : template <typename T>
112 2 : ValueCache<T>::ValueCache(const std::string & file_name,
113 : std::size_t in_dim,
114 : std::size_t max_leaf_size)
115 2 : : ValueCache(in_dim, max_leaf_size)
116 : {
117 2 : _persistent_storage_file = file_name;
118 :
119 : // if the persistent storage file exists and is readable, load it
120 2 : if (MooseUtils::checkFileReadable(*_persistent_storage_file,
121 : /*check_line_endings =*/false,
122 : /*throw_on_unreadable =*/false))
123 : {
124 1 : std::ifstream in_file(_persistent_storage_file->c_str());
125 1 : if (!in_file)
126 0 : mooseError("Failed to open '", *_persistent_storage_file, "' for reading.");
127 1 : dataLoad(in_file, *this, nullptr);
128 1 : }
129 2 : }
130 :
131 : template <typename T>
132 6 : ValueCache<T>::~ValueCache()
133 : {
134 : // if a persistent storage file was specified, write results back to it
135 6 : if (_persistent_storage_file.has_value())
136 : {
137 2 : std::ofstream out_file(_persistent_storage_file->c_str());
138 2 : if (!out_file)
139 0 : mooseWarning("Failed to open '", *_persistent_storage_file, "' for writing.");
140 2 : dataStore(out_file, *this, nullptr);
141 2 : }
142 6 : }
143 :
144 : template <typename T>
145 : void
146 120 : ValueCache<T>::insert(const std::vector<Real> & in_val, const T & out_val)
147 : {
148 : mooseAssert(in_val.size() == _in_dim, "Key dimensions do not match cache dimensions");
149 :
150 120 : auto id = size();
151 120 : _location_data.emplace_back(in_val, out_val);
152 :
153 : // update bounding box
154 120 : if (id == 0)
155 : {
156 : // first item is inserted
157 4 : _bbox.resize(_in_dim);
158 14 : for (const auto i : make_range(_in_dim))
159 10 : _bbox[i] = {in_val[i], in_val[i]};
160 : }
161 : else
162 246 : for (const auto i : make_range(_in_dim))
163 130 : _bbox[i] = {std::min(_bbox[i].first, in_val[i]), std::max(_bbox[i].second, in_val[i])};
164 :
165 : // do we have too many subtrees?
166 120 : if (_kd_tree && _kd_tree->getAllIndices().size() > _max_subtrees)
167 0 : _kd_tree = nullptr;
168 :
169 : // rebuild tree or add point
170 120 : if (!_kd_tree)
171 : {
172 4 : _kd_tree = std::make_unique<KdTreeT>(_in_dim, *this, _max_leaf_size);
173 : mooseAssert(_kd_tree != nullptr, "KDTree was not properly initialized.");
174 : }
175 : else
176 116 : _kd_tree->addPoints(id, id);
177 120 : }
178 :
179 : /**
180 : * Retrieve a single closest neighbor to in_val. Throws an exception if no values are stored.
181 : */
182 : template <typename T>
183 : std::tuple<const std::vector<Real> &, const T &, Real>
184 3 : ValueCache<T>::getNeighbor(const std::vector<Real> & in_val)
185 : {
186 : // throw an exception if this is called on an empty cache
187 3 : if (_location_data.empty())
188 1 : throw std::runtime_error("Attempting to retrieve a neighbor from an empty ValueCache.");
189 :
190 : // buffers for the kNN search
191 2 : nanoflann::KNNResultSet<Real> result_set(1);
192 : std::size_t return_index;
193 : Real distance;
194 :
195 : // kNN search
196 2 : result_set.init(&return_index, &distance);
197 2 : _kd_tree->findNeighbors(result_set, in_val.data());
198 :
199 2 : const auto & [location, data] = _location_data[return_index];
200 2 : return {std::cref(location), std::cref(data), distance};
201 : }
202 :
203 : /**
204 : * This function performs a search on the value cache and returns either the k-nearest neighbors or
205 : * the neighbors available if the cache size is less than k.
206 : */
207 : template <typename T>
208 : std::vector<std::tuple<const std::vector<Real> &, const T &, Real>>
209 4 : ValueCache<T>::getNeighbors(const std::vector<Real> & in_val, const std::size_t k)
210 : {
211 : // return early if no points are stored
212 4 : if (_location_data.empty())
213 1 : return {};
214 :
215 : // buffers for the kNN search
216 3 : nanoflann::KNNResultSet<Real> result_set(std::min(k, size()));
217 3 : std::vector<std::size_t> return_indices(std::min(k, size()));
218 3 : std::vector<Real> distances(std::min(k, size()));
219 :
220 : // kNN search
221 3 : result_set.init(return_indices.data(), distances.data());
222 3 : _kd_tree->findNeighbors(result_set, in_val.data());
223 :
224 : // prepare results to be returned
225 3 : std::vector<std::tuple<const std::vector<Real> &, const T &, Real>> nearest_neighbors;
226 10 : for (const auto i : index_range(result_set))
227 : {
228 7 : const auto & [location, data] = _location_data[return_indices[i]];
229 7 : nearest_neighbors.emplace_back(std::cref(location), std::cref(data), distances[i]);
230 : }
231 3 : return nearest_neighbors;
232 3 : }
233 :
234 : template <typename T>
235 : std::size_t
236 129 : ValueCache<T>::size()
237 : {
238 129 : return kdtree_get_point_count();
239 : }
240 :
241 : template <typename T>
242 : void
243 1 : ValueCache<T>::clear()
244 : {
245 1 : _location_data.clear();
246 1 : _kd_tree = nullptr;
247 1 : }
248 :
249 : template <typename T>
250 : void
251 1 : ValueCache<T>::rebuildTree()
252 : {
253 1 : if (_location_data.empty())
254 0 : return;
255 :
256 : // reset bounding box (must be done before the tree is built)
257 1 : _bbox.resize(_in_dim);
258 1 : const auto & location0 = _location_data[0].first;
259 4 : for (const auto i : make_range(_in_dim))
260 3 : _bbox[i] = {location0[i], location0[i]};
261 :
262 3 : for (const auto & pair : _location_data)
263 : {
264 2 : const auto & location = pair.first;
265 8 : for (const auto i : make_range(_in_dim))
266 6 : _bbox[i] = {std::min(_bbox[i].first, location[i]), std::max(_bbox[i].second, location[i])};
267 : }
268 :
269 : // build kd-tree
270 1 : _kd_tree = std::make_unique<KdTreeT>(_in_dim, *this, _max_leaf_size);
271 : mooseAssert(_kd_tree != nullptr, "KDTree was not properly initialized.");
272 : }
273 :
274 : template <typename T>
275 : std::size_t
276 134 : ValueCache<T>::kdtree_get_point_count() const
277 : {
278 134 : return _location_data.size();
279 : }
280 :
281 : template <typename T>
282 : Real
283 2231 : ValueCache<T>::kdtree_get_pt(const std::size_t idx, const std::size_t dim) const
284 : {
285 2231 : return _location_data[idx].first[dim];
286 : }
287 :
288 : template <typename T>
289 : template <class BBOX>
290 : bool
291 121 : ValueCache<T>::kdtree_get_bbox(BBOX & bb) const
292 : {
293 121 : if (_location_data.empty())
294 0 : return false;
295 :
296 : // return the bounding box incrementally built upon insertion
297 264 : for (const auto i : make_range(_in_dim))
298 143 : bb[i] = {_bbox[i].first, _bbox[i].second};
299 121 : return true;
300 : }
301 :
302 : template <typename T>
303 : inline void
304 2 : dataStore(std::ostream & stream, ValueCache<T> & c, void * context)
305 : {
306 2 : storeHelper(stream, c._location_data, context);
307 2 : }
308 :
309 : template <typename T>
310 : inline void
311 1 : dataLoad(std::istream & stream, ValueCache<T> & c, void * context)
312 : {
313 1 : loadHelper(stream, c._location_data, context);
314 1 : c.rebuildTree();
315 1 : }
|
