OGS 6.2.0-97-g4a610c866
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  {
37  OGS_FATAL(
38  "The d.o.f. table passed must be for one variable that has "
39  "only one component!");
40  }
41 }
42 
44  const unsigned num_components,
45  ExtrapolatableElementCollection const& extrapolatables,
46  const double t,
47  GlobalVector const& current_solution,
48  LocalToGlobalIndexMap const& dof_table)
49 {
50  auto const num_nodal_dof_result =
52 
53  std::vector<GlobalIndexType> ghost_indices;
54  { // Create num_components times version of ghost_indices arranged by
55  // location. For example for 3 components and ghost_indices {5,6,10} we
56  // compute {15, 16, 17, 18, 19, 20, 30, 31, 32}.
57  auto const& single_component_ghost_indices =
59  auto const single_component_ghost_indices_size =
60  single_component_ghost_indices.size();
61  ghost_indices.reserve(single_component_ghost_indices_size *
62  num_components);
63  for (unsigned i = 0; i < single_component_ghost_indices_size; ++i)
64  {
65  for (unsigned c = 0; c < num_components; ++c)
66  {
67  ghost_indices.push_back(
68  single_component_ghost_indices[i] * num_components + c);
69  }
70  }
71  }
72 
73  if (!_nodal_values ||
74 #ifdef USE_PETSC
75  _nodal_values->getLocalSize() + _nodal_values->getGhostSize()
76 #else
77  _nodal_values->size()
78 #endif
79  != static_cast<GlobalIndexType>(num_nodal_dof_result))
80  {
82  {num_nodal_dof_result, num_nodal_dof_result, &ghost_indices,
83  nullptr});
84  }
85  _nodal_values->setZero();
86 
87  // counts the writes to each nodal value, i.e., the summands in order to
88  // compute the average afterwards
89  auto counts =
91  counts->setZero();
92 
93  auto const size = extrapolatables.size();
94  for (std::size_t i = 0; i < size; ++i)
95  {
96  extrapolateElement(i, num_components, extrapolatables, t,
97  current_solution, dof_table, *counts);
98  }
100 
102  *counts);
103 }
104 
106  const unsigned num_components,
107  ExtrapolatableElementCollection const& extrapolatables,
108  const double t,
109  GlobalVector const& current_solution,
110  LocalToGlobalIndexMap const& dof_table)
111 {
112  auto const num_element_dof_result = static_cast<GlobalIndexType>(
113  _dof_table_single_component.size() * num_components);
114 
115  if (!_residuals || _residuals->size() != num_element_dof_result)
116  {
117 #ifndef USE_PETSC
118  _residuals.reset(new GlobalVector{num_element_dof_result});
119 #else
120  _residuals.reset(new GlobalVector{num_element_dof_result, false});
121 #endif
122  }
123 
124  if (static_cast<std::size_t>(num_element_dof_result) !=
125  extrapolatables.size() * num_components)
126  {
127  OGS_FATAL("mismatch in number of D.o.F.");
128  }
129 
130  auto const size = extrapolatables.size();
131  for (std::size_t i = 0; i < size; ++i)
132  {
133  calculateResidualElement(i, num_components, extrapolatables, t,
134  current_solution, dof_table);
135  }
137 }
138 
140  std::size_t const element_index,
141  const unsigned num_components,
142  ExtrapolatableElementCollection const& extrapolatables,
143  const double t,
144  GlobalVector const& current_solution,
145  LocalToGlobalIndexMap const& dof_table,
146  GlobalVector& counts)
147 {
148  auto const& integration_point_values =
149  extrapolatables.getIntegrationPointValues(
150  element_index, t, current_solution, dof_table,
152 
153  auto const& N_0 = extrapolatables.getShapeMatrix(element_index, 0);
154  auto const num_nodes = static_cast<unsigned>(N_0.cols());
155  auto const num_values =
156  static_cast<unsigned>(integration_point_values.size());
157 
158  if (num_values % num_components != 0)
159  {
160  OGS_FATAL(
161  "The number of computed integration point values is not divisable "
162  "by the number of num_components. Maybe the computed property is "
163  "not a %d-component vector for each integration point.",
164  num_components);
165  }
166 
167  // number of integration points in the element
168  const auto num_int_pts = num_values / num_components;
169 
170  if (num_int_pts < num_nodes)
171  {
172  OGS_FATAL(
173  "Least squares is not possible if there are more nodes than"
174  "integration points.");
175  }
176 
177  auto const pair_it_inserted = _qr_decomposition_cache.emplace(
178  std::make_pair(num_nodes, num_int_pts), CachedData{});
179 
180  auto& cached_data = pair_it_inserted.first->second;
181  if (pair_it_inserted.second)
182  {
183  DBUG("Computing new singular value decomposition");
184 
185  // interpolation_matrix * nodal_values = integration_point_values
186  // We are going to pseudo-invert this relation now using singular value
187  // decomposition.
188  auto& interpolation_matrix = cached_data.A;
189  interpolation_matrix.resize(num_int_pts, num_nodes);
190 
191  interpolation_matrix.row(0) = N_0;
192  for (unsigned int_pt = 1; int_pt < num_int_pts; ++int_pt)
193  {
194  auto const& shp_mat =
195  extrapolatables.getShapeMatrix(element_index, int_pt);
196  assert(shp_mat.cols() == num_nodes);
197 
198  // copy shape matrix to extrapolation matrix row-wise
199  interpolation_matrix.row(int_pt) = shp_mat;
200  }
201 
202  // JacobiSVD is extremely reliable, but fast only for small matrices.
203  // But we usually have small matrices and we don't compute very often.
204  // Cf.
205  // http://eigen.tuxfamily.org/dox/group__TopicLinearAlgebraDecompositions.html
206  //
207  // Decomposes interpolation_matrix = U S V^T.
208  Eigen::JacobiSVD<Eigen::MatrixXd> svd(
209  interpolation_matrix, Eigen::ComputeThinU | Eigen::ComputeThinV);
210 
211  auto const& S = svd.singularValues();
212  auto const& U = svd.matrixU();
213  auto const& V = svd.matrixV();
214 
215  // Compute and save the pseudo inverse V * S^{-1} * U^T.
216  auto const rank = svd.rank();
217  assert(rank == num_nodes);
218 
219  // cf. http://eigen.tuxfamily.org/dox/JacobiSVD_8h_source.html
220  cached_data.A_pinv.noalias() = V.leftCols(rank) *
221  S.head(rank).asDiagonal().inverse() *
222  U.leftCols(rank).transpose();
223  }
224  else if (cached_data.A.row(0) != N_0)
225  {
226  OGS_FATAL("The cached and the passed shapematrices differ.");
227  }
228 
229  auto const& global_indices =
230  _dof_table_single_component(element_index, 0).rows;
231 
232  if (num_components == 1)
233  {
234  auto const integration_point_values_vec =
235  MathLib::toVector(integration_point_values);
236 
237  // Apply the pre-computed pseudo-inverse.
238  Eigen::VectorXd const nodal_values =
239  cached_data.A_pinv * integration_point_values_vec;
240 
241  // TODO does that give rise to PETSc problems? E.g., writing to ghost
242  // nodes? Furthermore: Is ghost nodes communication necessary for PETSc?
243  _nodal_values->add(global_indices, nodal_values);
244  counts.add(global_indices,
245  std::vector<double>(global_indices.size(), 1.0));
246  }
247  else
248  {
249  auto const integration_point_values_mat = MathLib::toMatrix(
250  integration_point_values, num_components, num_int_pts);
251 
252  // Apply the pre-computed pseudo-inverse.
253  Eigen::MatrixXd const nodal_values =
254  cached_data.A_pinv * integration_point_values_mat.transpose();
255 
256  std::vector<GlobalIndexType> indices;
257  indices.reserve(num_components * global_indices.size());
258 
259  // _nodal_values is ordered location-wise
260  for (unsigned comp = 0; comp < num_components; ++comp)
261  {
262  for (auto i : global_indices)
263  {
264  indices.push_back(num_components * i + comp);
265  }
266  }
267 
268  // Nodal_values are passed as a raw pointer, because PETScVector and
269  // EigenVector implementations differ slightly.
270  _nodal_values->add(indices, nodal_values.data());
271  counts.add(indices, std::vector<double>(indices.size(), 1.0));
272  }
273 }
274 
276  std::size_t const element_index,
277  const unsigned num_components,
278  ExtrapolatableElementCollection const& extrapolatables,
279  const double t,
280  GlobalVector const& current_solution,
281  LocalToGlobalIndexMap const& dof_table)
282 {
283  auto const& int_pt_vals = extrapolatables.getIntegrationPointValues(
284  element_index, t, current_solution, dof_table,
286 
287  auto const num_values = static_cast<unsigned>(int_pt_vals.size());
288  if (num_values % num_components != 0)
289  {
290  OGS_FATAL(
291  "The number of computed integration point values is not divisable "
292  "by the number of num_components. Maybe the computed property is "
293  "not a %d-component vector for each integration point.",
294  num_components);
295  }
296 
297  // number of integration points in the element
298  const auto num_int_pts = num_values / num_components;
299 
300  const auto& global_indices =
301  _dof_table_single_component(element_index, 0).rows;
302  const auto num_nodes = static_cast<unsigned>(global_indices.size());
303 
304  auto const& interpolation_matrix =
305  _qr_decomposition_cache.find({num_nodes, num_int_pts})->second.A;
306 
307  Eigen::VectorXd nodal_vals_element(num_nodes);
308  auto const int_pt_vals_mat =
309  MathLib::toMatrix(int_pt_vals, num_components, num_int_pts);
310 
311  MathLib::LinAlg::setLocalAccessibleVector(
312  *_nodal_values); // For access in the for-loop.
313  for (unsigned comp = 0; comp < num_components; ++comp)
314  {
315  // filter nodal values of the current element
316  for (unsigned i = 0; i < num_nodes; ++i)
317  {
318  // TODO PETSc negative indices?
319  auto const idx = num_components * global_indices[i] + comp;
320  nodal_vals_element[i] = _nodal_values->get(idx);
321  }
322 
323  double const residual = (interpolation_matrix * nodal_vals_element -
324  int_pt_vals_mat.row(comp).transpose())
325  .squaredNorm();
326 
327  auto const eidx =
328  static_cast<GlobalIndexType>(num_components * element_index + comp);
329  // The residual is set to the root mean square value.
330  auto const root_mean_square = std::sqrt(residual / num_int_pts);
331  _residuals->set(eidx, root_mean_square);
332  }
333 }
334 
335 } // 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:63
GlobalMatrix::IndexType GlobalIndexType
std::vector< GlobalIndexType > const & getGhostIndices() const
Get ghost indices, forwarded from MeshComponentMap.
Stores a matrix and its Moore-Penrose pseudo-inverse.