OGS 6.1.0-1721-g6382411ad
LocalLinearLeastSquaresExtrapolator.cpp
Go to the documentation of this file.
1 
11 
12 #include <Eigen/SVD>
13 #include <logog/include/logog.hpp>
14 
16 #include "MathLib/LinAlg/LinAlg.h"
21 
22 namespace NumLib
23 {
25  NumLib::LocalToGlobalIndexMap const& dof_table)
26  : _dof_table_single_component(dof_table)
27 {
28  /* Note in case the following assertion fails:
29  * If you copied the extrapolation code, for your processes from
30  * somewhere, note that the code from the groundwater flow process might
31  * not suit your needs: It is a special case and is therefore most
32  * likely too simplistic. You better adapt the extrapolation code from
33  * some more advanced process, like the TES process.
34  */
35  if (dof_table.getNumberOfComponents() != 1)
36  OGS_FATAL(
37  "The d.o.f. table passed must be for one variable that has "
38  "only one component!");
39 }
40 
42  const unsigned num_components,
43  ExtrapolatableElementCollection const& extrapolatables,
44  const double t,
45  GlobalVector const& current_solution,
46  LocalToGlobalIndexMap const& dof_table)
47 {
48  auto const num_nodal_dof_result =
50 
51  std::vector<GlobalIndexType> ghost_indices;
52  { // Create num_components times version of ghost_indices arranged by
53  // location. For example for 3 components and ghost_indices {5,6,10} we
54  // compute {15, 16, 17, 18, 19, 20, 30, 31, 32}.
55  auto const& single_component_ghost_indices =
57  auto const single_component_ghost_indices_size =
58  single_component_ghost_indices.size();
59  ghost_indices.reserve(single_component_ghost_indices_size *
60  num_components);
61  for (unsigned i = 0; i < single_component_ghost_indices_size; ++i)
62  {
63  for (unsigned c = 0; c < num_components; ++c)
64  {
65  ghost_indices.push_back(
66  single_component_ghost_indices[i] * num_components + c);
67  }
68  }
69  }
70 
71  if (!_nodal_values ||
72 #ifdef USE_PETSC
73  _nodal_values->getLocalSize() + _nodal_values->getGhostSize()
74 #else
75  _nodal_values->size()
76 #endif
77  != static_cast<GlobalIndexType>(num_nodal_dof_result))
78  {
80  {num_nodal_dof_result, num_nodal_dof_result, &ghost_indices,
81  nullptr});
82  }
83  _nodal_values->setZero();
84 
85  // counts the writes to each nodal value, i.e., the summands in order to
86  // compute the average afterwards
87  auto counts =
89  counts->setZero();
90 
91  auto const size = extrapolatables.size();
92  for (std::size_t i = 0; i < size; ++i)
93  {
94  extrapolateElement(i, num_components, extrapolatables, t,
95  current_solution, dof_table, *counts);
96  }
98 
100  *counts);
101 }
102 
104  const unsigned num_components,
105  ExtrapolatableElementCollection const& extrapolatables,
106  const double t,
107  GlobalVector const& current_solution,
108  LocalToGlobalIndexMap const& dof_table)
109 {
110  auto const num_element_dof_result = static_cast<GlobalIndexType>(
111  _dof_table_single_component.size() * num_components);
112 
113  if (!_residuals || _residuals->size() != num_element_dof_result)
114  {
115 #ifndef USE_PETSC
116  _residuals.reset(new GlobalVector{num_element_dof_result});
117 #else
118  _residuals.reset(new GlobalVector{num_element_dof_result, false});
119 #endif
120  }
121 
122  if (static_cast<std::size_t>(num_element_dof_result) !=
123  extrapolatables.size() * num_components)
124  {
125  OGS_FATAL("mismatch in number of D.o.F.");
126  }
127 
128  auto const size = extrapolatables.size();
129  for (std::size_t i = 0; i < size; ++i)
130  {
131  calculateResidualElement(i, num_components, extrapolatables, t,
132  current_solution, dof_table);
133  }
135 }
136 
138  std::size_t const element_index,
139  const unsigned num_components,
140  ExtrapolatableElementCollection const& extrapolatables,
141  const double t,
142  GlobalVector const& current_solution,
143  LocalToGlobalIndexMap const& dof_table,
144  GlobalVector& counts)
145 {
146  auto const& integration_point_values =
147  extrapolatables.getIntegrationPointValues(
148  element_index, t, current_solution, dof_table,
150 
151  auto const& N_0 = extrapolatables.getShapeMatrix(element_index, 0);
152  auto const num_nodes = static_cast<unsigned>(N_0.cols());
153  auto const num_values =
154  static_cast<unsigned>(integration_point_values.size());
155 
156  if (num_values % num_components != 0)
157  OGS_FATAL(
158  "The number of computed integration point values is not divisable "
159  "by the number of num_components. Maybe the computed property is "
160  "not a %d-component vector for each integration point.",
161  num_components);
162 
163  // number of integration points in the element
164  const auto num_int_pts = num_values / num_components;
165 
166  if (num_int_pts < num_nodes)
167  OGS_FATAL(
168  "Least squares is not possible if there are more nodes than"
169  "integration points.");
170 
171  auto const pair_it_inserted = _qr_decomposition_cache.emplace(
172  std::make_pair(num_nodes, num_int_pts), CachedData{});
173 
174  auto& cached_data = pair_it_inserted.first->second;
175  if (pair_it_inserted.second)
176  {
177  DBUG("Computing new singular value decomposition");
178 
179  // interpolation_matrix * nodal_values = integration_point_values
180  // We are going to pseudo-invert this relation now using singular value
181  // decomposition.
182  auto& interpolation_matrix = cached_data.A;
183  interpolation_matrix.resize(num_int_pts, num_nodes);
184 
185  interpolation_matrix.row(0) = N_0;
186  for (unsigned int_pt = 1; int_pt < num_int_pts; ++int_pt)
187  {
188  auto const& shp_mat =
189  extrapolatables.getShapeMatrix(element_index, int_pt);
190  assert(shp_mat.cols() == num_nodes);
191 
192  // copy shape matrix to extrapolation matrix row-wise
193  interpolation_matrix.row(int_pt) = shp_mat;
194  }
195 
196  // JacobiSVD is extremely reliable, but fast only for small matrices.
197  // But we usually have small matrices and we don't compute very often.
198  // Cf.
199  // http://eigen.tuxfamily.org/dox/group__TopicLinearAlgebraDecompositions.html
200  //
201  // Decomposes interpolation_matrix = U S V^T.
202  Eigen::JacobiSVD<Eigen::MatrixXd> svd(
203  interpolation_matrix, Eigen::ComputeThinU | Eigen::ComputeThinV);
204 
205  auto const& S = svd.singularValues();
206  auto const& U = svd.matrixU();
207  auto const& V = svd.matrixV();
208 
209  // Compute and save the pseudo inverse V * S^{-1} * U^T.
210  auto const rank = svd.rank();
211  assert(rank == num_nodes);
212 
213  // cf. http://eigen.tuxfamily.org/dox/JacobiSVD_8h_source.html
214  cached_data.A_pinv.noalias() = V.leftCols(rank) *
215  S.head(rank).asDiagonal().inverse() *
216  U.leftCols(rank).transpose();
217  }
218  else if (cached_data.A.row(0) != N_0)
219  {
220  OGS_FATAL("The cached and the passed shapematrices differ.");
221  }
222 
223  auto const& global_indices =
224  _dof_table_single_component(element_index, 0).rows;
225 
226  if (num_components == 1)
227  {
228  auto const integration_point_values_vec =
229  MathLib::toVector(integration_point_values);
230 
231  // Apply the pre-computed pseudo-inverse.
232  Eigen::VectorXd const nodal_values =
233  cached_data.A_pinv * integration_point_values_vec;
234 
235  // TODO does that give rise to PETSc problems? E.g., writing to ghost
236  // nodes? Furthermore: Is ghost nodes communication necessary for PETSc?
237  _nodal_values->add(global_indices, nodal_values);
238  counts.add(global_indices,
239  std::vector<double>(global_indices.size(), 1.0));
240  }
241  else
242  {
243  auto const integration_point_values_mat = MathLib::toMatrix(
244  integration_point_values, num_components, num_int_pts);
245 
246  // Apply the pre-computed pseudo-inverse.
247  Eigen::MatrixXd const nodal_values =
248  cached_data.A_pinv * integration_point_values_mat.transpose();
249 
250  std::vector<GlobalIndexType> indices;
251  indices.reserve(num_components * global_indices.size());
252 
253  // _nodal_values is ordered location-wise
254  for (unsigned comp = 0; comp < num_components; ++comp)
255  {
256  for (auto i : global_indices)
257  {
258  indices.push_back(num_components * i + comp);
259  }
260  }
261 
262  // Nodal_values are passed as a raw pointer, because PETScVector and
263  // EigenVector implementations differ slightly.
264  _nodal_values->add(indices, nodal_values.data());
265  counts.add(indices, std::vector<double>(indices.size(), 1.0));
266  }
267 }
268 
270  std::size_t const element_index,
271  const unsigned num_components,
272  ExtrapolatableElementCollection const& extrapolatables,
273  const double t,
274  GlobalVector const& current_solution,
275  LocalToGlobalIndexMap const& dof_table)
276 {
277  auto const& int_pt_vals = extrapolatables.getIntegrationPointValues(
278  element_index, t, current_solution, dof_table,
280 
281  auto const num_values = static_cast<unsigned>(int_pt_vals.size());
282  if (num_values % num_components != 0)
283  OGS_FATAL(
284  "The number of computed integration point values is not divisable "
285  "by the number of num_components. Maybe the computed property is "
286  "not a %d-component vector for each integration point.",
287  num_components);
288 
289  // number of integration points in the element
290  const auto num_int_pts = num_values / num_components;
291 
292  const auto& global_indices =
293  _dof_table_single_component(element_index, 0).rows;
294  const auto num_nodes = static_cast<unsigned>(global_indices.size());
295 
296  auto const& interpolation_matrix =
297  _qr_decomposition_cache.find({num_nodes, num_int_pts})->second.A;
298 
299  Eigen::VectorXd nodal_vals_element(num_nodes);
300  auto const int_pt_vals_mat =
301  MathLib::toMatrix(int_pt_vals, num_components, num_int_pts);
302 
303  MathLib::LinAlg::setLocalAccessibleVector(
304  *_nodal_values); // For access in the for-loop.
305  for (unsigned comp = 0; comp < num_components; ++comp)
306  {
307  // filter nodal values of the current element
308  for (unsigned i = 0; i < num_nodes; ++i)
309  {
310  // TODO PETSc negative indices?
311  auto const idx = num_components * global_indices[i] + comp;
312  nodal_vals_element[i] = _nodal_values->get(idx);
313  }
314 
315  double const residual = (interpolation_matrix * nodal_vals_element -
316  int_pt_vals_mat.row(comp).transpose())
317  .squaredNorm();
318 
319  auto const eidx =
320  static_cast<GlobalIndexType>(num_components * element_index + comp);
321  // The residual is set to the root mean square value.
322  auto const root_mean_square = std::sqrt(residual / num_int_pts);
323  _residuals->set(eidx, root_mean_square);
324  }
325 }
326 
327 } // namespace NumLib
std::vector< double > _integration_point_values_cache
Avoids frequent reallocations.
void extrapolateElement(std::size_t const element_index, const unsigned num_components, ExtrapolatableElementCollection const &extrapolatables, const double t, GlobalVector const &current_solution, LocalToGlobalIndexMap const &dof_table, GlobalVector &counts)
Extrapolate one element.
void calculateResiduals(const unsigned num_components, ExtrapolatableElementCollection const &extrapolatables, const double t, GlobalVector const &current_solution, LocalToGlobalIndexMap const &dof_table) override
void finalizeAssembly(Matrix &)
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.
virtual std::vector< double > const & getIntegrationPointValues(std::size_t const id, const double t, GlobalVector const &current_solution, LocalToGlobalIndexMap const &dof_table, std::vector< double > &cache) const =0
void extrapolate(const unsigned num_components, ExtrapolatableElementCollection const &extrapolatables, const double t, GlobalVector const &current_solution, LocalToGlobalIndexMap const &dof_table) override
Extrapolates the given property from the given local assemblers.
virtual std::size_t size() const =0
Returns the number of elements whose properties shall be extrapolated.
void componentwiseDivide(MatrixOrVector &w, MatrixOrVector const &x, MatrixOrVector const &y)
Computes componentwise.
void calculateResidualElement(std::size_t const element_index, const unsigned num_components, ExtrapolatableElementCollection const &extrapolatables, const double t, GlobalVector const &current_solution, LocalToGlobalIndexMap const &dof_table)
Compute the residuals for one element.
NumLib::LocalToGlobalIndexMap const & _dof_table_single_component
DOF table used for writing to global vectors.
std::unique_ptr< GlobalVector > _residuals
extrapolation residuals
virtual Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(std::size_t const id, unsigned const integration_point) const =0
LocalLinearLeastSquaresExtrapolator(NumLib::LocalToGlobalIndexMap const &dof_table)
std::unique_ptr< GlobalVector > _nodal_values
extrapolated nodal values
std::map< std::pair< unsigned, unsigned >, CachedData > _qr_decomposition_cache
Eigen::Map< const Matrix > toMatrix(std::vector< double > const &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:71
#define OGS_FATAL(fmt,...)
Definition: Error.h:71
GlobalMatrix::IndexType GlobalIndexType
std::vector< GlobalIndexType > const & getGhostIndices() const
Get ghost indices, forwarded from MeshComponentMap.
Stores a matrix and its Moore-Penrose pseudo-inverse.