Nektar++
GlobalLinSysPETScStaticCond.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: GlobalLinSysPETScStaticCond.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: GlobalLinSys definition
32//
33///////////////////////////////////////////////////////////////////////////////
34
36
37#include <petscksp.h>
38#include <petscmat.h>
39#include <petscsys.h>
40
41using namespace std;
42
44{
45/**
46 * @class GlobalLinSysPETSc
47 *
48 * Solves a linear system using single- or multi-level static
49 * condensation.
50 */
51
52/**
53 * Registers the class with the Factory.
54 */
57 "PETScStaticCond", GlobalLinSysPETScStaticCond::create,
58 "PETSc static condensation.");
59
62 "PETScMultiLevelStaticCond", GlobalLinSysPETScStaticCond::create,
63 "PETSc multi-level static condensation.");
64
65/**
66 * For a matrix system of the form @f[
67 * \left[ \begin{array}{cc}
68 * \boldsymbol{A} & \boldsymbol{B}\\
69 * \boldsymbol{C} & \boldsymbol{D}
70 * \end{array} \right]
71 * \left[ \begin{array}{c} \boldsymbol{x_1}\\ \boldsymbol{x_2}
72 * \end{array}\right]
73 * = \left[ \begin{array}{c} \boldsymbol{y_1}\\ \boldsymbol{y_2}
74 * \end{array}\right],
75 * @f]
76 * where @f$\boldsymbol{D}@f$ and
77 * @f$(\boldsymbol{A-BD^{-1}C})@f$ are invertible, store and assemble
78 * a static condensation system, according to a given local to global
79 * mapping. #m_linSys is constructed by AssembleSchurComplement().
80 * @param mKey Associated matrix key.
81 * @param pLocMatSys LocalMatrixSystem
82 * @param locToGloMap Local to global mapping.
83 */
85 const GlobalLinSysKey &pKey, const std::weak_ptr<ExpList> &pExpList,
86 const std::shared_ptr<AssemblyMap> &pLocToGloMap)
87 : GlobalLinSys(pKey, pExpList, pLocToGloMap),
88 GlobalLinSysPETSc(pKey, pExpList, pLocToGloMap),
89 GlobalLinSysStaticCond(pKey, pExpList, pLocToGloMap)
90{
93 "This constructor is only valid when using static "
94 "condensation");
96 pLocToGloMap->GetGlobalSysSolnType(),
97 "The local to global map is not set up for the requested "
98 "solution type");
99}
100
101/**
102 *
103 */
105 const GlobalLinSysKey &pKey, const std::weak_ptr<ExpList> &pExpList,
106 const DNekScalBlkMatSharedPtr pSchurCompl,
108 const DNekScalBlkMatSharedPtr pInvD,
109 const std::shared_ptr<AssemblyMap> &pLocToGloMap,
110 const PreconditionerSharedPtr pPrecon)
111 : GlobalLinSys(pKey, pExpList, pLocToGloMap),
112 GlobalLinSysPETSc(pKey, pExpList, pLocToGloMap),
113 GlobalLinSysStaticCond(pKey, pExpList, pLocToGloMap)
114{
115 m_schurCompl = pSchurCompl;
116 m_BinvD = pBinvD;
117 m_C = pC;
118 m_invD = pInvD;
119 m_precon = pPrecon;
120}
121
122/**
123 *
124 */
126{
127}
128
130{
131 auto asmMap = m_locToGloMap.lock();
132
133 m_precon = CreatePrecon(asmMap);
134
135 // Allocate memory for top-level structure
136 SetupTopLevel(asmMap);
137
138 // Setup Block Matrix systems
139 int n, n_exp = m_expList.lock()->GetNumElmts();
140
141 // Build preconditioner
142 m_precon->BuildPreconditioner();
143
144 // Do transform of Schur complement matrix
145 int cnt = 0;
146 for (n = 0; n < n_exp; ++n)
147 {
149 {
150 DNekScalMatSharedPtr mat = m_schurCompl->GetBlock(n, n);
152 m_precon->TransformedSchurCompl(n, cnt, mat);
153 m_schurCompl->SetBlock(n, n, t);
154 cnt += mat->GetRows();
155 }
156 }
157
158 // Construct this level
159 Initialise(asmMap);
160}
161
162/**
163 * Assemble the schur complement matrix from the block matrices stored
164 * in #m_blkMatrices and the given local to global mapping information.
165 * @param locToGloMap Local to global mapping information.
166 */
168 AssemblyMapSharedPtr pLocToGloMap)
169{
170 int i, j, n, cnt, gid1, gid2, loc_lda;
171 NekDouble sign1, sign2, value;
172
173 const int nDirDofs = pLocToGloMap->GetNumGlobalDirBndCoeffs();
174
179 DNekScalMatSharedPtr loc_mat;
180
181 // Build precon again if we in multi-level static condensation (a
182 // bit of a hack)
184 {
186 m_precon->BuildPreconditioner();
187 }
188
189 // CALCULATE REORDERING MAPPING
190 CalculateReordering(pLocToGloMap->GetGlobalToUniversalBndMap(),
191 pLocToGloMap->GetGlobalToUniversalBndMapUnique(),
192 pLocToGloMap);
193
194 // SET UP VECTORS AND MATRIX
195 SetUpMatVec(pLocToGloMap->GetNumGlobalBndCoeffs(), nDirDofs);
196
197 // SET UP SCATTER OBJECTS
198 SetUpScatter();
199
200 // CONSTRUCT KSP OBJECT
202 m_expList.lock()->GetSession();
203 string variable = pLocToGloMap->GetVariable();
204 NekDouble IterativeSolverTolerance = NekConstants::kNekIterativeTol;
205 if (pSession->DefinesGlobalSysSolnInfo(variable,
206 "IterativeSolverTolerance"))
207 {
208 IterativeSolverTolerance = boost::lexical_cast<double>(
209 pSession->GetGlobalSysSolnInfo(variable, "IterativeSolverTolerance")
210 .c_str());
211 }
212 else if (pSession->DefinesParameter("IterativeSolverTolerance"))
213 {
214 IterativeSolverTolerance =
215 pSession->GetParameter("IterativeSolverTolerance");
216 }
217 SetUpSolver(IterativeSolverTolerance);
218
219 // If we are using the matrix multiplication shell don't try to
220 // populate the matrix.
222 {
223 return;
224 }
225
226 // POPULATE MATRIX
227 for (n = cnt = 0; n < m_schurCompl->GetNumberOfBlockRows(); ++n)
228 {
229 loc_mat = m_schurCompl->GetBlock(n, n);
230 loc_lda = loc_mat->GetRows();
231
232 for (i = 0; i < loc_lda; ++i)
233 {
234 gid1 = pLocToGloMap->GetLocalToGlobalBndMap(cnt + i) - nDirDofs;
235 sign1 = pLocToGloMap->GetLocalToGlobalBndSign(cnt + i);
236 if (gid1 >= 0)
237 {
238 int gid1ro = m_reorderedMap[gid1];
239 for (j = 0; j < loc_lda; ++j)
240 {
241 gid2 = pLocToGloMap->GetLocalToGlobalBndMap(cnt + j) -
242 nDirDofs;
243 sign2 = pLocToGloMap->GetLocalToGlobalBndSign(cnt + j);
244 if (gid2 >= 0)
245 {
246 int gid2ro = m_reorderedMap[gid2];
247 value = sign1 * sign2 * (*loc_mat)(i, j);
248 MatSetValue(m_matrix, gid1ro, gid2ro, value,
249 ADD_VALUES);
250 }
251 }
252 }
253 }
254 cnt += loc_lda;
255 }
256
257 // ASSEMBLE MATRIX
258 MatAssemblyBegin(m_matrix, MAT_FINAL_ASSEMBLY);
259 MatAssemblyEnd(m_matrix, MAT_FINAL_ASSEMBLY);
260}
261
263 unsigned int n)
264{
265 DNekScalBlkMatSharedPtr schurComplBlock;
266 DNekScalMatSharedPtr localMat = m_schurCompl->GetBlock(n, n);
267 unsigned int nbdry = localMat->GetRows();
268 unsigned int nblks = 1;
269 unsigned int esize[1] = {nbdry};
270
272 nblks, nblks, esize, esize);
273 schurComplBlock->SetBlock(0, 0, localMat);
274
275 return schurComplBlock;
276}
277
279 int scLevel, [[maybe_unused]] Array<OneD, NekDouble> &F_bnd)
280{
281 if (scLevel == 0)
282 {
283 // When matrices are supplied to the constructor at the top
284 // level, the preconditioner is never set up.
285 if (!m_precon)
286 {
288 m_precon->BuildPreconditioner();
289 }
290 }
291}
292
293/**
294 * @brief Solve linear system using PETSc.
295 *
296 * The general strategy being a PETSc solve is to:
297 *
298 * - Copy values into the PETSc vector #m_b
299 * - Solve the system #m_ksp and place result into #m_x.
300 * - Scatter results back into #m_locVec using #m_ctx scatter object.
301 * - Copy from #m_locVec to output array #pOutput.
302 */
304 const int pNumRows, const Array<OneD, const NekDouble> &pInput,
305 Array<OneD, NekDouble> &pOutput, const AssemblyMapSharedPtr &locToGloMap,
306 const int pNumDir)
307{
308 const int nHomDofs = pNumRows - pNumDir;
309
311 {
312 m_precon = CreatePrecon(locToGloMap);
313 m_precon->BuildPreconditioner();
314 }
315
316 Array<OneD, NekDouble> Glo(pNumRows);
317 locToGloMap->AssembleBnd(pInput, Glo);
318
319 // Populate RHS vector from input
320 VecSetValues(m_b, nHomDofs, &m_reorderedMap[0], &Glo[pNumDir],
321 INSERT_VALUES);
322
323 // Assemble RHS vector
324 VecAssemblyBegin(m_b);
325 VecAssemblyEnd(m_b);
326
327 // Do system solve
328 KSPSolve(m_ksp, m_b, m_x);
329
330 KSPConvergedReason reason;
331 KSPGetConvergedReason(m_ksp, &reason);
332 ASSERTL0(reason > 0, "PETSc solver diverged, reason is: " +
333 std::string(KSPConvergedReasons[reason]));
334
335 // Scatter results to local vector
336 VecScatterBegin(m_ctx, m_x, m_locVec, INSERT_VALUES, SCATTER_FORWARD);
337 VecScatterEnd(m_ctx, m_x, m_locVec, INSERT_VALUES, SCATTER_FORWARD);
338
339 // Copy results into output vector
340 PetscScalar *tmp;
341 VecGetArray(m_locVec, &tmp);
342 Vmath::Vcopy(nHomDofs, tmp, 1, &Glo[pNumDir], 1);
343 Vmath::Zero(pNumDir, Glo, 1);
344 locToGloMap->GlobalToLocalBnd(Glo, pOutput);
345 VecRestoreArray(m_locVec, &tmp);
346}
347
350{
351 m_precon->DoTransformBasisToLowEnergy(pInOut);
352}
353
356{
357 m_precon->DoTransformCoeffsFromLowEnergy(pInOut);
358}
359
361 const Array<OneD, NekDouble> &pInput, Array<OneD, NekDouble> &pOutput)
362{
363 m_precon->DoTransformCoeffsToLowEnergy(pInput, pOutput);
364}
365
366/**
367 * @brief Apply matrix-vector multiplication using local approach and
368 * the assembly map.
369 *
370 * @param input Vector input.
371 * @param output Result of multiplication.
372 *
373 * @todo This can possibly be made faster by using the sparse
374 * block-matrix multiplication code from the iterative elastic
375 * systems.
376 */
379{
380 auto asmMap = m_locToGloMap.lock();
381
382 int nLocBndDofs = asmMap->GetNumLocalBndCoeffs();
383 int nDirDofs = asmMap->GetNumGlobalDirBndCoeffs();
384
385 NekVector<NekDouble> in(nLocBndDofs), out(nLocBndDofs);
386 asmMap->GlobalToLocalBnd(input, in.GetPtr(), nDirDofs);
387 out = (*m_schurCompl) * in;
388 asmMap->AssembleBnd(out.GetPtr(), output, nDirDofs);
389}
390
392 const GlobalLinSysKey &mkey, const std::weak_ptr<ExpList> &pExpList,
393 const DNekScalBlkMatSharedPtr pSchurCompl,
395 const DNekScalBlkMatSharedPtr pInvD,
396 const std::shared_ptr<AssemblyMap> &l2gMap)
397{
400 mkey, pExpList, pSchurCompl, pBinvD, pC, pInvD, l2gMap, m_precon);
401 sys->Initialise(l2gMap);
402 return sys;
403}
404} // namespace Nektar::MultiRegions
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
A global linear system.
Definition: GlobalLinSys.h:70
const std::weak_ptr< ExpList > m_expList
Local Matrix System.
Definition: GlobalLinSys.h:122
const GlobalLinSysKey m_linSysKey
Key associated with this linear system.
Definition: GlobalLinSys.h:120
PreconditionerSharedPtr CreatePrecon(AssemblyMapSharedPtr asmMap)
Create a preconditioner object from the parameters defined in the supplied assembly map.
void Initialise(const std::shared_ptr< AssemblyMap > &pLocToGloMap)
Definition: GlobalLinSys.h:203
GlobalSysSolnType GetGlobalSysSolnType() const
Return the associated solution type.
A PETSc global linear system.
Vec m_x
PETSc vector objects used for local storage.
PETScMatMult m_matMult
Enumerator to select matrix multiplication type.
std::vector< int > m_reorderedMap
Reordering that takes universal IDs to a unique row in the PETSc matrix.
void SetUpScatter()
Set up PETSc local (equivalent to Nektar++ global) and global (equivalent to universal) scatter maps.
void SetUpSolver(NekDouble tolerance)
Set up KSP solver object.
void SetUpMatVec(int nGlobal, int nDir)
Construct PETSc matrix and vector handles.
VecScatter m_ctx
PETSc scatter context that takes us between Nektar++ global ordering and PETSc vector ordering.
void CalculateReordering(const Array< OneD, const int > &glo2uniMap, const Array< OneD, const int > &glo2unique, const AssemblyMapSharedPtr &pLocToGloMap)
Calculate a reordering of universal IDs for PETSc.
KSP m_ksp
KSP object that represents solver system.
void v_CoeffsFwdTransform(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput) override
void v_AssembleSchurComplement(std::shared_ptr< AssemblyMap > locToGloMap) override
Assemble the Schur complement matrix.
static GlobalLinSysSharedPtr create(const GlobalLinSysKey &pLinSysKey, const std::weak_ptr< ExpList > &pExpList, const std::shared_ptr< AssemblyMap > &pLocToGloMap)
Creates an instance of this class.
void v_SolveLinearSystem(const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir) override
Solve linear system using PETSc.
GlobalLinSysStaticCondSharedPtr v_Recurse(const GlobalLinSysKey &mkey, const std::weak_ptr< ExpList > &pExpList, const DNekScalBlkMatSharedPtr pSchurCompl, const DNekScalBlkMatSharedPtr pBinvD, const DNekScalBlkMatSharedPtr pC, const DNekScalBlkMatSharedPtr pInvD, const std::shared_ptr< AssemblyMap > &locToGloMap) override
void v_CoeffsBwdTransform(Array< OneD, NekDouble > &pInOut) override
void v_DoMatrixMultiply(const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output) override
Apply matrix-vector multiplication using local approach and the assembly map.
DNekScalBlkMatSharedPtr v_GetStaticCondBlock(unsigned int n) override
Retrieves a the static condensation block matrices from n-th expansion using the matrix key provided ...
void v_PreSolve(int scLevel, Array< OneD, NekDouble > &F_bBnd) override
GlobalLinSysPETScStaticCond(const GlobalLinSysKey &mkey, const std::weak_ptr< ExpList > &pExpList, const std::shared_ptr< AssemblyMap > &locToGloMap)
Constructor for full direct matrix solve.
void v_BasisFwdTransform(Array< OneD, NekDouble > &pInOut) override
DNekScalBlkMatSharedPtr m_schurCompl
Block Schur complement matrix.
std::weak_ptr< AssemblyMap > m_locToGloMap
Local to global map.
void SetupTopLevel(const std::shared_ptr< AssemblyMap > &locToGloMap)
Set up the storage for the Schur complement or the top level of the multi-level Schur complement.
DNekScalBlkMatSharedPtr m_BinvD
Block matrix.
DNekScalBlkMatSharedPtr m_C
Block matrix.
DNekScalBlkMatSharedPtr m_invD
Block matrix.
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
Array< OneD, DataType > & GetPtr()
Definition: NekVector.cpp:217
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< GlobalLinSysStaticCond > GlobalLinSysStaticCondSharedPtr
GlobalLinSysFactory & GetGlobalLinSysFactory()
std::shared_ptr< Preconditioner > PreconditionerSharedPtr
Definition: GlobalLinSys.h:58
std::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
Definition: AssemblyMap.h:50
std::shared_ptr< GlobalLinSysPETScStaticCond > GlobalLinSysPETScStaticCondSharedPtr
static const NekDouble kNekIterativeTol
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:79
double NekDouble
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
STL namespace.