chombo-discharge
Loading...
Searching...
No Matches
CD_LeastSquaresImplem.H
1/* chombo-discharge
2 * Copyright © 2021 SINTEF Energy Research.
3 * Please refer to Copyright.txt and LICENSE in the chombo-discharge root directory.
4 */
5
12#ifndef CD_LeastSquaresImplem_H
13#define CD_LeastSquaresImplem_H
14
15// Our includes
16#include <CD_LaPackUtils.H>
17#include <CD_LeastSquares.H>
18#include <CD_MultiIndex.H>
19#include <CD_NamespaceHeader.H>
20
21inline Vector<Real>
23{
24 Vector<Real> ret(a_displacements.size(), 1.0);
25
26 if (a_pow > 0) {
27 for (int i = 0; i < a_displacements.size(); i++) {
28 const RealVect& d = a_displacements[i];
29 const Real len = d.vectorLength();
30 const Real w = std::pow(len, -a_pow);
31
32 ret[i] = w;
33 }
34 }
35
36 return ret;
37}
38
39template <typename T>
63
64template <typename T>
74 const int a_order)
75{
76 // TLDR: This routine does a two-level least squares solve of a system of equations, each equation describing
77 // a Taylor series expansion up to specified order (a_order). This routine "solves" this overdetermined system
78 // using weighted least squares so that we can obtain terms in the Taylor series. The user will have specified
79 // the desired unknowns (a_derivs) to be returned in stencil form as well as any terms that were already
80 // known (a_knownTerms).
81 //
82 // The least squares solve is done by LaPack, and so the least squares system matrix is filled in Fortran
83 // order. The actual "solve" consists of using the Moore-Penrose pseudoinverse and once we have inverted
84 // the system we can form stencils for each of the specified terms in the Taylor series.
85 //
86 // This routine does not contain much code, but can be difficult to understand.
87
88 CH_assert(a_order > 0);
91 CH_assert(a_fineVofs.size() == a_fineWeights.size());
92 CH_assert(a_coarVofs.size() == a_coarWeights.size());
93
94 // Initialize return stuff
96 for (IVSIterator ivsIt(a_derivs); ivsIt.ok(); ++ivsIt) {
98 }
99
100 if (a_derivs.numPts() > 0) {
101
102 // This is because some unknowns (rows) can be eliminated.
104 const int Kfine = a_fineDisplacements.size();
105 const int Kcoar = a_coarDisplacements.size();
106 const int K = Kfine + Kcoar;
107
109
110 if (K < M) {
111 MayDay::Abort("LeastSquares::computeDualLevelStencils -- not enough equations to achieve desired order!");
112 }
113 if (!isect.isEmpty()) {
115 "LeastSquares::computeDualLevelStencils - you have specified the same terms as both unknown and known");
116 }
117
118 // Build the A-matrix in column major order (this is what Fortran wants) so we can use LaPackUtils::computePseudoInverse.
119 // ----------------------------------------------------------------------------------------------------------------------
120 // If we have an (unweighted) full system then our system A*x = b is
121 //
122 // A x b
123 // === === ===
124 // [1 (x-x0) (x-x0)^2 ...] [f(x) ] [f(x0)]
125 // [1 (x-x1) (x-x1)^2 ...] [df/dx ] [f(x1)]
126 // [1 (x-x2) (x-x2)^2 ...] [d^2f/dx^2] = [f(x2)]
127 // [: : : ] [ : ] [ : ]
128 // [: : : ] [ : ] [ : ]
129 //
130 // Extensions to 2D/3D simply use multi-index notation, and weights are simply multiplied into each row, i.e. we
131 // solve (w*A) * x = (w*b). Inverting the system gives x = [w * A^+ * w] * b where A^+ is the Moore-Penrose
132 // pseudoinverse. We put the result of [w * A^+ * w] into a stencil.
133 //
134 // Note that columns can be eliminated through knownTerms, in which case we remove unknowns (i.e., rows) from the system.
135 // This will also correspond to a modification of the right-hand side, but the required modifications are not accesible
136 // in this routine, and so the user will have to make sense of them.
137
138 int i = 0; // Exists just because we fill memory linearly.
139 Vector<T> linA(K * M, 0.0); // Equal to (w*A)
140 Vector<T> linAplus(M * K, 0.0); // Equal to (w*A)^+
141
142 // Loop over column
143 for (MultiIndex mi(a_order); mi.ok(); ++mi) {
144 if (!a_knownTerms.contains(mi.getCurrentIndex())) {
145
146 // Fill column, write the system using the fine vofs first, then the coarse vofs.
147 for (int k = 0; k < K; k++) {
148 if (k < Kfine) {
149 linA[i] = a_fineWeights[k] * mi.pow(a_fineDisplacements[k]) / mi.factorial();
150 }
151 else {
152 linA[i] = a_coarWeights[k - Kfine] * mi.pow(a_coarDisplacements[k - Kfine]) / mi.factorial();
153 }
154 i++;
155 }
156 }
157 }
158
159 // Compute the pseudo-inverse.
160 const bool foundSVD = LaPackUtils::computePseudoInverse(linAplus.stdVector(), linA.stdVector(), K, M);
161
162 if (foundSVD) {
163 // When we have eliminated rows in the linear system we can't use MultiIndex to map directly, this occurs because
164 // if we eliminated unknowns (rows) our system can be something like (if we eliminated term f(x)):
165 //
166 // [(x-x0) (x-x0)^2 ...] [df/dx ] [f(x0) - f(x)]
167 // [(x-x1) (x-x1)^2 ...] [d^2f/dx^2] [f(x1) - f(x)]
168 // [(x-x2) (x-x2)^2 ...] [ ] = [f(x2) - f(x)]
169 // [ : : ] [ : ] [ : ]
170 // [ : : ] [ : ] [ : ]
171 //
172 // The MultiIndex won't know about this, and will (correctly!), believe that the first row corresponds to
173 // multi-index (0,0,0). Since that is truly not the case, we map the rows in A to multi-indices and use that to
174 // identify the row in (w*A)^+ that corresponds to a specific unknown in the Taylor series. This is what happens below.
176 int row = 0;
177 for (MultiIndex mi(a_order); mi.ok(); ++mi) {
178
179 // This is the order in which A was built.
180 if (!a_knownTerms.contains(mi.getCurrentIndex())) {
181 rowMap.emplace(mi.getCurrentIndex(), row);
182 row++;
183 }
184 }
185
186 // Recall that linAplus is M*K so the stride is always M, starting at some specified row.
187 for (IVSIterator ivsIt(a_derivs); ivsIt.ok(); ++ivsIt) {
188 const IntVect deriv = ivsIt();
189
190 if (rowMap.find(ivsIt()) != rowMap.end()) {
191 row = rowMap.at(ivsIt());
192 }
193 else {
194 MayDay::Error("LeastSquares::computeDualLevelStencils -- map is out of range but this shouldn't happen!");
195 }
196
198
199 sten.first.clear();
200 sten.second.clear();
201
202 // Map the pseudoinverse into something that is usable by a stencil. Note that linAplus is (w*A)^+, but we want
203 // the term [(w*A)^+ * w], so we also need to multiply in the weights here (because the right-hand side was also weighted).
204 for (int k = 0; k < K; k++) {
205 const int idx = row + k * M;
206 if (k < Kfine) {
207 sten.first.add(a_fineVofs[k], a_fineWeights[k] * linAplus[idx]);
208 }
209 else {
210 sten.second.add(a_coarVofs[k - Kfine], a_coarWeights[k - Kfine] * linAplus[idx]);
211 }
212 }
213 }
214 }
215 else {
216 MayDay::Warning("LeastSquares::computeDualLevelStencils - could not perform singular value decomposition");
217 }
218 }
219
220 return ret;
221}
222
223#include <CD_NamespaceFooter.H>
224
225#endif
Interface to some LaPack routines.
Agglomeration of various routines that use least squares computations.
static std::map< IntVect, std::pair< VoFStencil, VoFStencil > > computeDualLevelStencils(const IntVectSet &a_derivs, const IntVectSet &a_knownTerms, const Vector< VolIndex > &a_fineVofs, const Vector< VolIndex > &a_coarVofs, const Vector< RealVect > &a_fineDisplacements, const Vector< RealVect > &a_coarDisplacements, const int a_p, const int a_order)
Compute a least squares interpolation to a specified order. This version separates the stencils into ...
Definition CD_LeastSquaresImplem.H:41
static int getTaylorExpansionSize(const int a_order)
Get the size of a Taylor expansion for a given order.
Definition CD_LeastSquares.cpp:416
static Vector< Real > makeDiagWeights(const Vector< RealVect > &a_displacements, const int a_pow)
Create a list of weights. This routine returns a list of diagonal weights for a least squares system....
Definition CD_LeastSquaresImplem.H:22
SpaceDim multi-index type for use with higher order Taylor series.
Definition CD_MultiIndex.H:30
Base class for a tracer particle solver. This solver can advance particles in a pre-defined velocity ...
Definition CD_TracerParticleSolver.H:37
TracerParticleSolver()
Default constructor.
Definition CD_TracerParticleSolverImplem.H:25
bool computePseudoInverse(std::vector< double > &a_linAplus, const std::vector< double > &a_linA, const int &a_M, const int &a_N)
Compute the pseudoinverse of matrix through singular value decomposition.
Definition CD_LaPackUtils.cpp:163