OGS
NumLib::LocalLinearLeastSquaresExtrapolator Class Reference

Detailed Description

Extrapolator doing a linear least squares extrapolation locally in each element.

The results of the element-wise least squares are averaged at each node.

Note

The described procedure is not least-squares optimal on the whole mesh! But the residuals are computed from the actual interpolation result that is being returned by this class.

The number of integration points in each element must be greater than or equal to the number of nodes of that element. This restriction is due to the use of the least squares which requires an exact or overdetermined equation system.

Definition at line 38 of file LocalLinearLeastSquaresExtrapolator.h.

#include <LocalLinearLeastSquaresExtrapolator.h>

Inheritance diagram for NumLib::LocalLinearLeastSquaresExtrapolator:
[legend]
Collaboration diagram for NumLib::LocalLinearLeastSquaresExtrapolator:
[legend]

Classes

struct  CachedData
 Stores a matrix and its Moore-Penrose pseudo-inverse. More...
 

Public Member Functions

 LocalLinearLeastSquaresExtrapolator (NumLib::LocalToGlobalIndexMap const &dof_table)
 
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.
 
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
 
GlobalVector const & getNodalValues () const override
 
GlobalVector const & getElementResiduals () const override
 
- Public Member Functions inherited from NumLib::Extrapolator
virtual ~Extrapolator ()=default
 

Private Member Functions

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.
 
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.
 

Private Attributes

std::unique_ptr< GlobalVector_nodal_values
 extrapolated nodal values
 
std::unique_ptr< GlobalVector_residuals
 extrapolation residuals
 
NumLib::LocalToGlobalIndexMap const & _dof_table_single_component
 DOF table used for writing to global vectors.
 
std::vector< double > _integration_point_values_cache
 Avoids frequent reallocations.
 
std::map< std::pair< unsigned, unsigned >, CachedData_qr_decomposition_cache
 

Constructor & Destructor Documentation

◆ LocalLinearLeastSquaresExtrapolator()

NumLib::LocalLinearLeastSquaresExtrapolator::LocalLinearLeastSquaresExtrapolator ( NumLib::LocalToGlobalIndexMap const & dof_table)
explicit

Constructs a new instance.

Note
The dof_table must point to a d.o.f. table for one single-component variable.

Definition at line 25 of file LocalLinearLeastSquaresExtrapolator.cpp.

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 {
39 "The d.o.f. table passed must be for one variable that has only "
40 "one component!");
41 }
42}
#define OGS_FATAL(...)
Definition Error.h:26
NumLib::LocalToGlobalIndexMap const & _dof_table_single_component
DOF table used for writing to global vectors.

References NumLib::LocalToGlobalIndexMap::getNumberOfGlobalComponents(), and OGS_FATAL.

Member Function Documentation

◆ calculateResidualElement()

void NumLib::LocalLinearLeastSquaresExtrapolator::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 )
private

Compute the residuals for one element.

Definition at line 281 of file LocalLinearLeastSquaresExtrapolator.cpp.

288{
289 auto const& int_pt_vals = extrapolatables.getIntegrationPointValues(
290 element_index, t, x, dof_table, _integration_point_values_cache);
291
292 auto const num_values = static_cast<unsigned>(int_pt_vals.size());
293 if (num_values % num_components != 0)
294 {
295 OGS_FATAL(
296 "The number of computed integration point values is not divisible "
297 "by the number of num_components. Maybe the computed property is "
298 "not a {:d}-component vector for each integration point.",
299 num_components);
300 }
301
302 // number of integration points in the element
303 const auto num_int_pts = num_values / num_components;
304
305 const auto& global_indices =
306 _dof_table_single_component(element_index, 0).rows;
307 const auto num_nodes = static_cast<unsigned>(global_indices.size());
308
309 auto const& interpolation_matrix =
310 _qr_decomposition_cache.find({num_nodes, num_int_pts})->second.A;
311
312 Eigen::VectorXd nodal_vals_element(num_nodes);
313 auto const int_pt_vals_mat =
314 MathLib::toMatrix(int_pt_vals, num_components, num_int_pts);
315
317 *_nodal_values); // For access in the for-loop.
318 for (unsigned comp = 0; comp < num_components; ++comp)
319 {
320 // filter nodal values of the current element
321 for (unsigned i = 0; i < num_nodes; ++i)
322 {
323 // TODO PETSc negative indices?
324 auto const idx = num_components * global_indices[i] + comp;
325 nodal_vals_element[i] = _nodal_values->get(idx);
326 }
327
328 double const residual = (interpolation_matrix * nodal_vals_element -
329 int_pt_vals_mat.row(comp).transpose())
330 .squaredNorm();
331
332 auto const eidx =
333 static_cast<GlobalIndexType>(num_components * element_index + comp);
334 // The residual is set to the root mean square value.
335 auto const root_mean_square = std::sqrt(residual / num_int_pts);
336 _residuals->set(eidx, root_mean_square);
337 }
338}
GlobalMatrix::IndexType GlobalIndexType
std::unique_ptr< GlobalVector > _nodal_values
extrapolated nodal values
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
void setLocalAccessibleVector(PETScVector const &x)
Definition LinAlg.cpp:27
Eigen::Map< const Matrix > toMatrix(std::vector< double > const &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)

References _dof_table_single_component, _integration_point_values_cache, _nodal_values, _qr_decomposition_cache, _residuals, NumLib::ExtrapolatableElementCollection::getIntegrationPointValues(), OGS_FATAL, MathLib::LinAlg::setLocalAccessibleVector(), and MathLib::toMatrix().

Referenced by calculateResiduals().

◆ calculateResiduals()

void NumLib::LocalLinearLeastSquaresExtrapolator::calculateResiduals ( const unsigned num_components,
ExtrapolatableElementCollection const & extrapolatables,
const double t,
std::vector< GlobalVector * > const & x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const & dof_table )
overridevirtual

Computes residuals from the extrapolation of the given property.

The residuals are computed as element values.

Precondition
extrapolate() must have been called before with the same arguments.

The computed residuals are root-mean-square of the difference between the integration point values obtained from the local assemblers and the extrapolation results when interpolated back to the integration points again.

Implements NumLib::Extrapolator.

Definition at line 107 of file LocalLinearLeastSquaresExtrapolator.cpp.

113{
114 auto const num_element_dof_result = static_cast<GlobalIndexType>(
115 _dof_table_single_component.size() * num_components);
116
117 if (!_residuals || _residuals->size() != num_element_dof_result)
118 {
119#ifndef USE_PETSC
120 _residuals.reset(new GlobalVector{num_element_dof_result});
121#else
122 _residuals.reset(new GlobalVector{num_element_dof_result, false});
123#endif
124 }
125
126 if (static_cast<std::size_t>(num_element_dof_result) !=
127 extrapolatables.size() * num_components)
128 {
129 OGS_FATAL("mismatch in number of D.o.F.");
130 }
131
132 auto const size = extrapolatables.size();
133 for (std::size_t i = 0; i < size; ++i)
134 {
135 calculateResidualElement(i, num_components, extrapolatables, t, x,
136 dof_table);
137 }
139}
Global vector based on Eigen vector.
Definition EigenVector.h:25
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 finalizeAssembly(PETScMatrix &A)
Definition LinAlg.cpp:170
constexpr int size(int const displacement_dim)
Vectorized tensor size for given displacement dimension.

References _dof_table_single_component, _residuals, calculateResidualElement(), MathLib::LinAlg::finalizeAssembly(), OGS_FATAL, NumLib::LocalToGlobalIndexMap::size(), and NumLib::ExtrapolatableElementCollection::size().

◆ extrapolate()

void NumLib::LocalLinearLeastSquaresExtrapolator::extrapolate ( const unsigned num_components,
ExtrapolatableElementCollection const & extrapolatables,
const double t,
std::vector< GlobalVector * > const & x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const & dof_table )
overridevirtual

Extrapolates the given property from the given local assemblers.

Implements NumLib::Extrapolator.

Definition at line 44 of file LocalLinearLeastSquaresExtrapolator.cpp.

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 }
102
104 *counts);
105}
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.
std::vector< GlobalIndexType > const & getGhostIndices() const
Get ghost indices, forwarded from MeshComponentMap.
void componentwiseDivide(PETScVector &w, PETScVector const &x, PETScVector const &y)
Definition LinAlg.cpp:81

References _dof_table_single_component, _nodal_values, MathLib::LinAlg::componentwiseDivide(), NumLib::LocalToGlobalIndexMap::dofSizeWithoutGhosts(), extrapolateElement(), MathLib::LinAlg::finalizeAssembly(), NumLib::LocalToGlobalIndexMap::getGhostIndices(), and NumLib::ExtrapolatableElementCollection::size().

◆ extrapolateElement()

void NumLib::LocalLinearLeastSquaresExtrapolator::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 )
private

Extrapolate one element.

Definition at line 141 of file LocalLinearLeastSquaresExtrapolator.cpp.

149{
150 auto const& integration_point_values =
151 extrapolatables.getIntegrationPointValues(
152 element_index, t, x, dof_table, _integration_point_values_cache);
153
154 // Empty vector means to ignore the values and not to change the counts.
155 if (integration_point_values.empty())
156 {
157 return;
158 }
159
160 auto const& N_0 = extrapolatables.getShapeMatrix(element_index, 0);
161 auto const num_nodes = static_cast<unsigned>(N_0.cols());
162 auto const num_values =
163 static_cast<unsigned>(integration_point_values.size());
164
165 if (num_values % num_components != 0)
166 {
167 OGS_FATAL(
168 "The number of computed integration point values is not divisible "
169 "by the number of num_components. Maybe the computed property is "
170 "not a {:d}-component vector for each integration point.",
171 num_components);
172 }
173
174 // number of integration points in the element
175 const auto num_int_pts = num_values / num_components;
176
177 if (num_int_pts < num_nodes)
178 {
179 OGS_FATAL(
180 "Least squares is not possible if there are more nodes than "
181 "integration points.");
182 }
183
184 auto const pair_it_inserted = _qr_decomposition_cache.emplace(
185 std::make_pair(num_nodes, num_int_pts), CachedData{});
186
187 auto& cached_data = pair_it_inserted.first->second;
188 if (pair_it_inserted.second)
189 {
190 DBUG("Computing new singular value decomposition");
191
192 // interpolation_matrix * nodal_values = integration_point_values
193 // We are going to pseudo-invert this relation now using singular value
194 // decomposition.
195 auto& interpolation_matrix = cached_data.A;
196 interpolation_matrix.resize(num_int_pts, num_nodes);
197
198 interpolation_matrix.row(0) = N_0;
199 for (unsigned int_pt = 1; int_pt < num_int_pts; ++int_pt)
200 {
201 auto const& shp_mat =
202 extrapolatables.getShapeMatrix(element_index, int_pt);
203 assert(shp_mat.cols() == num_nodes);
204
205 // copy shape matrix to extrapolation matrix row-wise
206 interpolation_matrix.row(int_pt) = shp_mat;
207 }
208
209 // JacobiSVD is extremely reliable, but fast only for small matrices.
210 // But we usually have small matrices and we don't compute very often.
211 // Cf.
212 // http://eigen.tuxfamily.org/dox/group__TopicLinearAlgebraDecompositions.html
213 //
214 // Decomposes interpolation_matrix = U S V^T.
215 Eigen::JacobiSVD<Eigen::MatrixXd> svd(
216 interpolation_matrix, Eigen::ComputeThinU | Eigen::ComputeThinV);
217
218 auto const& S = svd.singularValues();
219 auto const& U = svd.matrixU();
220 auto const& V = svd.matrixV();
221
222 // Compute and save the pseudo inverse V * S^{-1} * U^T.
223 auto const rank = svd.rank();
224 assert(rank == num_nodes);
225
226 // cf. http://eigen.tuxfamily.org/dox/JacobiSVD_8h_source.html
227 cached_data.A_pinv.noalias() = V.leftCols(rank) *
228 S.head(rank).asDiagonal().inverse() *
229 U.leftCols(rank).transpose();
230 }
231 else if (cached_data.A.row(0) != N_0)
232 {
233 OGS_FATAL("The cached and the passed shapematrices differ.");
234 }
235
236 auto const& global_indices =
237 _dof_table_single_component(element_index, 0).rows;
238
239 if (num_components == 1)
240 {
241 auto const integration_point_values_vec =
242 MathLib::toVector(integration_point_values);
243
244 // Apply the pre-computed pseudo-inverse.
245 Eigen::VectorXd const nodal_values =
246 cached_data.A_pinv * integration_point_values_vec;
247
248 // TODO does that give rise to PETSc problems? E.g., writing to ghost
249 // nodes? Furthermore: Is ghost nodes communication necessary for PETSc?
250 _nodal_values->add(global_indices, nodal_values);
251 counts.add(global_indices,
252 std::vector<double>(global_indices.size(), 1.0));
253 }
254 else
255 {
256 auto const integration_point_values_mat = MathLib::toMatrix(
257 integration_point_values, num_components, num_int_pts);
258
259 // Apply the pre-computed pseudo-inverse.
260 Eigen::MatrixXd const nodal_values =
261 cached_data.A_pinv * integration_point_values_mat.transpose();
262
263 std::vector<GlobalIndexType> indices;
264 indices.reserve(num_components * global_indices.size());
265
266 // _nodal_values is ordered location-wise
267 for (unsigned comp = 0; comp < num_components; ++comp)
268 {
269 transform(cbegin(global_indices), cend(global_indices),
270 back_inserter(indices),
271 [&](auto const i) { return num_components * i + comp; });
272 }
273
274 // Nodal_values are passed as a raw pointer, because PETScVector and
275 // EigenVector implementations differ slightly.
276 _nodal_values->add(indices, nodal_values.data());
277 counts.add(indices, std::vector<double>(indices.size(), 1.0));
278 }
279}
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
void add(IndexType rowId, double v)
add entry
Definition EigenVector.h:76
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.

References _dof_table_single_component, _integration_point_values_cache, _nodal_values, _qr_decomposition_cache, MathLib::EigenVector::add(), DBUG(), NumLib::ExtrapolatableElementCollection::getIntegrationPointValues(), NumLib::ExtrapolatableElementCollection::getShapeMatrix(), OGS_FATAL, MathLib::toMatrix(), and MathLib::toVector().

Referenced by extrapolate().

◆ getElementResiduals()

GlobalVector const & NumLib::LocalLinearLeastSquaresExtrapolator::getElementResiduals ( ) const
inlineoverridevirtual

Returns the extrapolation residuals.

Todo
Maybe write directly to a MeshProperty.

Implements NumLib::Extrapolator.

Definition at line 77 of file LocalLinearLeastSquaresExtrapolator.h.

78 {
79 return *_residuals;
80 }

References _residuals.

◆ getNodalValues()

GlobalVector const & NumLib::LocalLinearLeastSquaresExtrapolator::getNodalValues ( ) const
inlineoverridevirtual

Returns the extrapolated nodal values.

Todo
Maybe write directly to a MeshProperty.

Implements NumLib::Extrapolator.

Definition at line 72 of file LocalLinearLeastSquaresExtrapolator.h.

73 {
74 return *_nodal_values;
75 }

References _nodal_values.

Member Data Documentation

◆ _dof_table_single_component

NumLib::LocalToGlobalIndexMap const& NumLib::LocalLinearLeastSquaresExtrapolator::_dof_table_single_component
private

DOF table used for writing to global vectors.

Definition at line 104 of file LocalLinearLeastSquaresExtrapolator.h.

Referenced by calculateResidualElement(), calculateResiduals(), extrapolate(), and extrapolateElement().

◆ _integration_point_values_cache

std::vector<double> NumLib::LocalLinearLeastSquaresExtrapolator::_integration_point_values_cache
private

Avoids frequent reallocations.

Definition at line 107 of file LocalLinearLeastSquaresExtrapolator.h.

Referenced by calculateResidualElement(), and extrapolateElement().

◆ _nodal_values

std::unique_ptr<GlobalVector> NumLib::LocalLinearLeastSquaresExtrapolator::_nodal_values
private

extrapolated nodal values

Definition at line 100 of file LocalLinearLeastSquaresExtrapolator.h.

Referenced by calculateResidualElement(), extrapolate(), extrapolateElement(), and getNodalValues().

◆ _qr_decomposition_cache

std::map<std::pair<unsigned, unsigned>, CachedData> NumLib::LocalLinearLeastSquaresExtrapolator::_qr_decomposition_cache
private

Maps (#nodes, #int_pts) to (N_0, QR decomposition), where N_0 is the shape matrix of the first integration point.

Note
It is assumed that the pair (#nodes, #int_pts) uniquely identifies the set of all shape matrices N for a mesh element (i.e., only N, not dN/dx etc.).
Todo
Add the element dimension as identifying criterion, or change to typeid.

Definition at line 129 of file LocalLinearLeastSquaresExtrapolator.h.

Referenced by calculateResidualElement(), and extrapolateElement().

◆ _residuals

std::unique_ptr<GlobalVector> NumLib::LocalLinearLeastSquaresExtrapolator::_residuals
private

extrapolation residuals

Definition at line 101 of file LocalLinearLeastSquaresExtrapolator.h.

Referenced by calculateResidualElement(), calculateResiduals(), and getElementResiduals().


The documentation for this class was generated from the following files: