G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsINSPreconditioners.hpp
Go to the documentation of this file.
1
12#pragma once
13#include <gsIncompressibleFlow/src/gsINSPreconditioners.h>
14
15namespace gismo
16{
17
18template <class T, int MatOrder>
19typename gsINSPreconditioner<T, MatOrder>::uPtr gsINSPreconditioner<T, MatOrder>::make(std::string precType, const std::map<std::string, gsSparseMatrix<T, MatOrder> >& mat, const gsOptionList& opt)
20{
21 if (precType == "LSC_FdiagEqual")
22 return gsINSBlockPrecondLSC<T, MatOrder, gsINSPrecondBlockF<T, MatOrder> >::make(mat, opt);
23 else if (precType == "LSC_Fdiag")
24 return gsINSBlockPrecondLSC<T, MatOrder, gsINSPrecondBlockFdiag<T, MatOrder> >::make(mat, opt);
25 else if (precType == "LSC_Fwhole")
26 return gsINSBlockPrecondLSC<T, MatOrder, gsINSPrecondBlockFwhole<T, MatOrder> >::make(mat, opt);
27 else if (precType == "LSC_Fmod")
28 return gsINSBlockPrecondLSC<T, MatOrder, gsINSPrecondBlockFmod<T, MatOrder> >::make(mat, opt);
29 else if (precType == "PCD_FdiagEqual")
30 return gsINSBlockPrecondPCD<T, MatOrder, gsINSPrecondBlockF<T, MatOrder> >::make(mat, opt);
31 else if (precType == "PCD_Fdiag")
32 return gsINSBlockPrecondPCD<T, MatOrder, gsINSPrecondBlockFdiag<T, MatOrder> >::make(mat, opt);
33 else if (precType == "PCD_Fwhole")
34 return gsINSBlockPrecondPCD<T, MatOrder, gsINSPrecondBlockFwhole<T, MatOrder> >::make(mat, opt);
35 else if (precType == "PCD_Fmod")
36 return gsINSBlockPrecondPCD<T, MatOrder, gsINSPrecondBlockFmod<T, MatOrder> >::make(mat, opt);
37 else if (precType == "PCDmod_FdiagEqual")
38 return gsINSBlockPrecondPCDmod<T, MatOrder, gsINSPrecondBlockF<T, MatOrder> >::make(mat, opt);
39 else if (precType == "PCDmod_Fdiag")
40 return gsINSBlockPrecondPCDmod<T, MatOrder, gsINSPrecondBlockFdiag<T, MatOrder> >::make(mat, opt);
41 else if (precType == "PCDmod_Fwhole")
42 return gsINSBlockPrecondPCDmod<T, MatOrder, gsINSPrecondBlockFwhole<T, MatOrder> >::make(mat, opt);
43 else if (precType == "PCDmod_Fmod")
44 return gsINSBlockPrecondPCDmod<T, MatOrder, gsINSPrecondBlockFmod<T, MatOrder> >::make(mat, opt);
45 else if (precType == "AL_Fwhole")
46 return gsINSBlockPrecondAL<T, MatOrder, gsINSPrecondBlockFwhole<T, MatOrder> >::make(mat, opt);
47 else if (precType == "AL_Fmod")
48 return gsINSBlockPrecondAL<T, MatOrder, gsINSPrecondBlockFmod<T, MatOrder> >::make(mat, opt);
49 else if (precType == "SIMPLE_FdiagEqual")
50 return gsINSBlockPrecondSIMPLE<T, MatOrder, gsINSPrecondBlockF<T, MatOrder> >::make(mat, opt);
51 else if (precType == "SIMPLE_Fdiag")
52 return gsINSBlockPrecondSIMPLE<T, MatOrder, gsINSPrecondBlockFdiag<T, MatOrder> >::make(mat, opt);
53 else if (precType == "SIMPLE_Fwhole")
54 return gsINSBlockPrecondSIMPLE<T, MatOrder, gsINSPrecondBlockFwhole<T, MatOrder> >::make(mat, opt);
55 else if (precType == "SIMPLE_Fmod")
56 return gsINSBlockPrecondSIMPLE<T, MatOrder, gsINSPrecondBlockFmod<T, MatOrder> >::make(mat, opt);
57 else if (precType == "SIMPLER_FdiagEqual")
58 return gsINSBlockPrecondSIMPLER<T, MatOrder, gsINSPrecondBlockF<T, MatOrder> >::make(mat, opt);
59 else if (precType == "SIMPLER_Fdiag")
60 return gsINSBlockPrecondSIMPLER<T, MatOrder, gsINSPrecondBlockFdiag<T, MatOrder> >::make(mat, opt);
61 else if (precType == "SIMPLER_Fwhole")
62 return gsINSBlockPrecondSIMPLER<T, MatOrder, gsINSPrecondBlockFwhole<T, MatOrder> >::make(mat, opt);
63 else if (precType == "SIMPLER_Fmod")
64 return gsINSBlockPrecondSIMPLER<T, MatOrder, gsINSPrecondBlockFmod<T, MatOrder> >::make(mat, opt);
65 else if (precType == "MSIMPLER_FdiagEqual")
66 return gsINSBlockPrecondMSIMPLER<T, MatOrder, gsINSPrecondBlockF<T, MatOrder> >::make(mat, opt);
67 else if (precType == "MSIMPLER_Fdiag")
68 return gsINSBlockPrecondMSIMPLER<T, MatOrder, gsINSPrecondBlockFdiag<T, MatOrder> >::make(mat, opt);
69 else if (precType == "MSIMPLER_Fwhole")
70 return gsINSBlockPrecondMSIMPLER<T, MatOrder, gsINSPrecondBlockFwhole<T, MatOrder> >::make(mat, opt);
71 else if (precType == "MSIMPLER_Fmod")
72 return gsINSBlockPrecondMSIMPLER<T, MatOrder, gsINSPrecondBlockFmod<T, MatOrder> >::make(mat, opt);
73 else if (precType == "StokesDiag_FdiagEqual")
74 return gsBlockPrecondStokes<T, MatOrder, gsINSPrecondBlockF<T, MatOrder> >::make(mat, opt);
75 else if (precType == "StokesTriang_FdiagEqual")
76 return gsBlockPrecondStokesTriang<T, MatOrder, gsINSPrecondBlockF<T, MatOrder> >::make(mat, opt);
77 /*else if (precType == "SchurEx_LSC_FdiagEqual")
78 return uwbINSBlockPrecondSchurEx<T, gsINSPrecondBlockF<T>, gsINSPrecondSchurLSC<T> >::make(mat, opt);
79 else if (precType == "SchurEx_SIMPLE_FdiagEqual")
80 return uwbINSBlockPrecondSchurEx<T, gsINSPrecondBlockF<T>, uwbINSPrecondSchurSIMPLE<T> >::make(mat, opt);
81 else if (precType == "SchurEx_MSIMPLER_FdiagEqual")
82 return uwbINSBlockPrecondSchurEx<T, gsINSPrecondBlockF<T>, gsINSPrecondSchurMSIMPLER<T> >::make(mat, opt);*/
83 else
84 {
85 gsInfo << "Invalid preconditioner type, using LSC_FdiagEqual.\n";
86 return gsINSBlockPrecondLSC<T, MatOrder, gsINSPrecondBlockF<T, MatOrder> >::make(mat, opt);
87 }
88}
89
90
91template <class T, int MatOrder>
92gsOptionList gsINSPreconditioner<T, MatOrder>::defaultOptions()
93{
94 gsOptionList opt;
95
96 opt.addInt("dim", "Problem dimension", 0);
97 opt.addReal("visc", "Viscosity", 0.0);
98 opt.addInt("udofs", "Number of velocity dofs", 0);
99 opt.addInt("pdofs", "Number of pressure dofs", 0);
100
101 opt.addInt("pcd_bcType", "Type of BCs for Fp, Ap in PCD preconditioner", 3);
102 opt.addReal("gamma", "Parameter gamma for AL", 1.0);
103 opt.addReal("alphaP", "Pressure relaxation parameter for SIMPLE-type", 1.0);
104 opt.addSwitch("lumpingM", "Use lumped diagonal mass matrices in preconditioners", true);
105 opt.addSwitch("lumpingA", "Use lumped diagonal blockA in SIMPLE-type preconditioners", false);
106 opt.addSwitch("pcd_assembAp", "Assemble pressure Poisson matrix for definition of Ap in PCD", true);
107 opt.addSwitch("pcd_assembFp", "Assemble pressure Poisson matrix for definition of Fp in PCD", true);
108
109 opt.addSwitch("iter", "Subsystems iteratively", false);
110 opt.addInt("maxIt", "Max number of iterations for inner solves", 10);
111 opt.addReal("tol", "Stopping tolerance for inner solves", 1e-2);
112 opt.addInt("fill", "Fill factor for ILUT precond (inner solves)", 2);
113 opt.addReal("dropTol", "Drop tolerance for ILUT precond (inner solves)", 1e-8);
114
115 return opt;
116}
117
118
119template <class T, int MatOrder>
120void gsINSBlockPrecondBase<T, MatOrder>::apply(const gsMatrix<T> & input, gsMatrix<T> & x) const
121{
122 int uSize = m_Finv->rows();
123 int pSize = m_Sinv->rows();
124
125 x.resize(uSize + pSize, 1);
126
127 gsMatrix<T> tmpX;
128 m_Sinv->apply(input.bottomRows(pSize), tmpX);
129 x.bottomRows(pSize) = tmpX;
130
131 gsMatrix<T> tmp;
132 m_Bt->apply(x.bottomRows(pSize), tmp);
133 tmp = input.topRows(uSize) - tmp;
134
135 m_Finv->apply(tmp, tmpX);
136 x.topRows(uSize) = tmpX;
137}
138
139
140template <class T, int MatOrder, class BlockFType>
141void gsINSBlockPrecondAL<T, MatOrder, BlockFType>::fillALgammaPart_into(gsSparseMatrix<T, MatOrder>& matGammaPart, gsMatrix<T>& rhsGammaPart, const std::map<std::string, gsSparseMatrix<T, MatOrder> >& mat, const gsMatrix<T>& rhs, const gsOptionList& opt)
142{
143 gsInfo << "Filling the extra part of Agamma matrix... ";
144
145 int dim = opt.getInt("dim");
146 int udofs = opt.getInt("udofs");
147 int pdofs = opt.getInt("pdofs");
148 int uSize = dim * udofs;
149 int numDofs = uSize + pdofs;
150 real_t gamma = opt.getReal("gamma");
151
152 const gsSparseMatrix<T, MatOrder>& matNS = mat.at("matNS");
153 const gsSparseMatrix<T, MatOrder>& presM = mat.at("matMp");
154
155 gsSparseMatrix<T, MatOrder> presMinv(pdofs, pdofs); // approximation of pressure mass matrix inverse
156 presMinv.setIdentity();
157 for (int i = 0; i < pdofs; i++)
158 presMinv.coeffRef(i, i) = 1 / presM.coeff(i, i);
159
160 gsSparseMatrix<T, MatOrder> Mgamma = gamma * matNS.block(0, uSize, uSize, pdofs) * presMinv * matNS.block(uSize, 0, pdofs, uSize);
161
162 gsVector<index_t> nonZerosVector;
163 nonZerosVector.setZero(numDofs);
164 nonZerosVector.topRows(uSize) = getNnzVectorPerOuter(Mgamma);
165
166 matGammaPart.resize(numDofs, numDofs);
167 matGammaPart.reserve(nonZerosVector);
168
169 for (int outer = 0; outer < Mgamma.outerSize(); ++outer)
170 for (typename gsSparseMatrix<T, MatOrder>::InnerIterator it(Mgamma, outer); it; ++it)
171 matGammaPart.insert(it.row(), it.col()) = it.value();
172
173 matGammaPart.makeCompressed();
174
175 rhsGammaPart.setZero(numDofs, 1);
176 rhsGammaPart.topRows(uSize) = gamma * matNS.block(0, uSize, uSize, pdofs) * presMinv * rhs.bottomRows(pdofs);
177
178 gsInfo << "Done.\n";
179}
180
181template <class T, int MatOrder, class BlockFType>
182void gsINSBlockPrecondAL<T, MatOrder, BlockFType>::fillALmodifSystem_into(gsSparseMatrix<T, MatOrder>& matGamma, gsMatrix<T>& rhsGamma, const std::map<std::string, gsSparseMatrix<T, MatOrder> >& mat, const gsMatrix<T>& rhs, const gsOptionList& opt)
183{
184 gsInfo << "Filling the Agamma matrix... ";
185
186 gsSparseMatrix<T, MatOrder> matGammaPart;
187 gsMatrix<T> rhsGammaPart;
188
189 fillALgammaPart_into(matGammaPart, rhsGammaPart, mat, rhs, opt);
190
191 matGamma = mat.at("matNS") + matGammaPart;
192 rhsGamma = rhs + rhsGammaPart;
193
194 gsInfo << "Done.\n";
195}
196
197
198template <class T, int MatOrder, class BlockFType>
199void gsINSBlockPrecondSIMPLE<T, MatOrder, BlockFType>::apply(const gsMatrix<T> & input, gsMatrix<T> & x) const
200{
201 int uSize = m_Finv->rows();
202 int pSize = m_Sinv->rows();
203
204 x.resize(uSize + pSize, 1);
205
206 gsMatrix<T> uStar;
207 m_Finv->apply(input.topRows(uSize), uStar);
208
209 gsMatrix<T> tmp1;
210 tmp1.noalias() = m_B * uStar;
211
212 gsMatrix<T> dp;
213 m_Sinv->apply(input.bottomRows(pSize) - tmp1, dp);
214
215 gsMatrix<T> tmp2;
216 m_Bt->apply(dp, tmp2);
217
218 x.topRows(uSize).noalias() = uStar - m_Dinv * tmp2;
219 x.bottomRows(pSize) = m_alphaP * dp;
220}
221
222
223template <class T, int MatOrder, class BlockFType>
224void gsINSBlockPrecondSIMPLER<T, MatOrder, BlockFType>::apply(const gsMatrix<T> & input, gsMatrix<T> & x) const
225{
226 int uSize = m_Finv->rows();
227 int pSize = m_Sinv->rows();
228
229 x.resize(uSize + pSize, 1);
230
231 // pressure estimate
232 gsMatrix<T> tmp1;
233 tmp1.noalias() = input.bottomRows(pSize) - m_BDinv * input.topRows(uSize);
234 gsMatrix<T> pStar;
235 m_Sinv->apply(tmp1, pStar);
236
237 // velocity solve
238 gsMatrix<T> BtpStar, uStar;
239 m_Bt->apply(pStar, BtpStar);
240 m_Finv->apply(input.topRows(uSize) - BtpStar, uStar);
241
242 // pressure correction
243 gsMatrix<T> dp;
244 m_Sinv->apply(input.bottomRows(pSize) - m_B * uStar, dp);
245
246 gsMatrix<T> tmp2;
247 m_Bt->apply(dp, tmp2);
248
249 x.topRows(uSize).noalias() = uStar - m_Dinv * tmp2;
250 x.bottomRows(pSize) = pStar + m_alphaP * dp;
251}
252
253
254template <class T, int MatOrder, class BlockFType>
255void gsINSBlockPrecondMSIMPLER<T, MatOrder, BlockFType>::apply(const gsMatrix<T> & input, gsMatrix<T> & x) const
256{
257 int uSize = m_Finv->rows();
258 int pSize = m_Sinv->rows();
259
260 x.resize(uSize + pSize, 1);
261
262 // pressure estimate
263 gsMatrix<T> tmp1;
264 tmp1.noalias() = input.bottomRows(pSize) - m_BMinv * input.topRows(uSize);
265 gsMatrix<T> pStar;
266 m_Sinv->apply(tmp1, pStar);
267
268 // velocity solve
269 gsMatrix<T> BtpStar, uStar;
270 m_Bt->apply(pStar, BtpStar);
271 m_Finv->apply(input.topRows(uSize) - BtpStar, uStar);
272
273 // pressure correction
274 gsMatrix<T> dp;
275 m_Sinv->apply(input.bottomRows(pSize) - m_B * uStar, dp);
276
277 gsMatrix<T> tmp2;
278 m_Bt->apply(dp, tmp2);
279
280 x.topRows(uSize).noalias() = uStar - m_velMinv * tmp2;
281 x.bottomRows(pSize) = pStar + m_alphaP * dp;
282}
283
284
285template <class T, int MatOrder, class BlockFType>
286void gsBlockPrecondStokes<T, MatOrder, BlockFType>::apply(const gsMatrix<T> & input, gsMatrix<T> & x) const
287{
288 int uSize = m_Finv->rows();
289 int pSize = m_Sinv->rows();
290
291 x.resize(uSize + pSize, 1);
292
293 gsMatrix<T> tmpX;
294 m_Sinv->apply(input.bottomRows(pSize), tmpX);
295 x.bottomRows(pSize) = tmpX;
296
297 m_Finv->apply(input.topRows(uSize), tmpX);
298 x.topRows(uSize) = tmpX;
299}
300
301// ===============================================================================
302
303
304
305} // namespace gismo
gismo::gsBlockPrecondStokesTriang
Block triangular preconditioner for the Stokes problem.
Definition gsINSPreconditioners.h:782
gismo::gsBlockPrecondStokes
Block diagonal preconditioner for the Stokes problem.
Definition gsINSPreconditioners.h:704
gismo::gsBlockPrecondStokes::apply
virtual void apply(const gsMatrix< T > &input, gsMatrix< T > &x) const
Apply the preconditioner. Computes the vector \fx = P^{-1} y\f.
Definition gsINSPreconditioners.hpp:286
gismo::gsINSBlockPrecondAL
Augmented Lagrangian preconditioner.
Definition gsINSPreconditioners.h:377
gismo::gsINSBlockPrecondAL::fillALmodifSystem_into
static void fillALmodifSystem_into(gsSparseMatrix< T, MatOrder > &matGamma, gsMatrix< T > &rhsGamma, const std::map< std::string, gsSparseMatrix< T, MatOrder > > &mat, const gsMatrix< T > &rhs, const gsOptionList &opt)
Fill the Augmented Lagrangian linear system.
Definition gsINSPreconditioners.hpp:182
gismo::gsINSBlockPrecondAL::fillALgammaPart_into
static void fillALgammaPart_into(gsSparseMatrix< T, MatOrder > &matGammaPart, gsMatrix< T > &rhsGammaPart, const std::map< std::string, gsSparseMatrix< T, MatOrder > > &mat, const gsMatrix< T > &rhs, const gsOptionList &opt)
Fill the extra part of the Augmented Lagrangian linear system.
Definition gsINSPreconditioners.hpp:141
gismo::gsINSBlockPrecondBase::apply
virtual void apply(const gsMatrix< T > &input, gsMatrix< T > &x) const
Apply the preconditioner. Computes the vector \fx = P^{-1} y\f.
Definition gsINSPreconditioners.hpp:120
gismo::gsINSBlockPrecondLSC
Least-squares commutator preconditioner.
Definition gsINSPreconditioners.h:191
gismo::gsINSBlockPrecondMSIMPLER
MSIMPLER preconditioner.
Definition gsINSPreconditioners.h:613
gismo::gsINSBlockPrecondMSIMPLER::apply
virtual void apply(const gsMatrix< T > &input, gsMatrix< T > &x) const
Apply the preconditioner. Computes the vector \fx = P^{-1} y\f.
Definition gsINSPreconditioners.hpp:255
gismo::gsINSBlockPrecondPCD
Pressure convection-diffusion preconditioner.
Definition gsINSPreconditioners.h:253
gismo::gsINSBlockPrecondPCDmod
Modified pressure convection-diffusion preconditioner.
Definition gsINSPreconditioners.h:315
gismo::gsINSBlockPrecondSIMPLER
SIMPLER preconditioner.
Definition gsINSPreconditioners.h:544
gismo::gsINSBlockPrecondSIMPLER::apply
virtual void apply(const gsMatrix< T > &input, gsMatrix< T > &x) const
Apply the preconditioner. Computes the vector \fx = P^{-1} y\f.
Definition gsINSPreconditioners.hpp:224
gismo::gsINSBlockPrecondSIMPLE
SIMPLE preconditioner.
Definition gsINSPreconditioners.h:457
gismo::gsINSBlockPrecondSIMPLE::apply
virtual void apply(const gsMatrix< T > &input, gsMatrix< T > &x) const
Apply the preconditioner. Computes the vector \fx = P^{-1} y\f.
Definition gsINSPreconditioners.hpp:199
gismo::gsMatrix
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition gsMatrix.h:41
gismo::gsOptionList
Class which holds a list of parameters/options, and provides easy access to them.
Definition gsOptionList.h:33
gismo::gsOptionList::addInt
void addInt(const std::string &label, const std::string &desc, const index_t &value)
Adds a option named label, with description desc and value value.
Definition gsOptionList.cpp:201
gismo::gsOptionList::getInt
const index_t & getInt(const std::string &label) const
Reads value for option label from options.
Definition gsOptionList.cpp:37
gismo::gsOptionList::getReal
Real getReal(const std::string &label) const
Reads value for option label from options.
Definition gsOptionList.cpp:44
gismo::gsOptionList::addReal
void addReal(const std::string &label, const std::string &desc, const Real &value)
Adds a option named label, with description desc and value value.
Definition gsOptionList.cpp:211
gismo::gsOptionList::addSwitch
void addSwitch(const std::string &label, const std::string &desc, const bool &value)
Adds a option named label, with description desc and value value.
Definition gsOptionList.cpp:235
gismo::gsSparseMatrix
Sparse matrix class, based on gsEigen::SparseMatrix.
Definition gsSparseMatrix.h:139
gismo::gsVector
A vector with arbitrary coefficient type and fixed or dynamic size.
Definition gsVector.h:37
gsInfo
#define gsInfo
Definition gsDebug.h:43
gismo
The G+Smo namespace, containing all definitions for the library.
gismo::getNnzVectorPerOuter
gsVector< index_t > getNnzVectorPerOuter(const gsSparseMatrix< T, MatOrder > &mat)
Get a vector of nonzero entries per outer index (row or column depending on the matrix storage order)...
Definition gsFlowUtils.h:116