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