Teko Version of the Day
Loading...
Searching...
No Matches
Teko_SIMPLEPreconditionerFactory.cpp
1/*
2// @HEADER
3//
4// ***********************************************************************
5//
6// Teko: A package for block and physics based preconditioning
7// Copyright 2010 Sandia Corporation
8//
9// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
10// the U.S. Government retains certain rights in this software.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric C. Cyr (eccyr@sandia.gov)
40//
41// ***********************************************************************
42//
43// @HEADER
44
45*/
46
47#include "Teko_SIMPLEPreconditionerFactory.hpp"
48
49#include "Teko_Utilities.hpp"
50#include "Teko_InverseFactory.hpp"
51#include "Teko_BlockLowerTriInverseOp.hpp"
52#include "Teko_BlockUpperTriInverseOp.hpp"
53#include <stdexcept>
54#ifdef TEKO_HAVE_EPETRA
55#include "Teko_DiagonalPreconditionerFactory.hpp"
56#endif
57
58#include "Teuchos_Time.hpp"
59
60using Teuchos::RCP;
61
62namespace Teko {
63namespace NS {
64
65// Constructor definition
66SIMPLEPreconditionerFactory
67 ::SIMPLEPreconditionerFactory(const RCP<InverseFactory> & inverse,
68 double alpha)
69 : invVelFactory_(inverse), invPrsFactory_(inverse), alpha_(alpha), fInverseType_(Diagonal), useMass_(false)
70{ }
71
72SIMPLEPreconditionerFactory
73 ::SIMPLEPreconditionerFactory(const RCP<InverseFactory> & invVFact,
74 const RCP<InverseFactory> & invPFact,
75 double alpha)
76 : invVelFactory_(invVFact), invPrsFactory_(invPFact), alpha_(alpha), fInverseType_(Diagonal), useMass_(false)
77{ }
78
79SIMPLEPreconditionerFactory::SIMPLEPreconditionerFactory()
80 : alpha_(1.0), fInverseType_(Diagonal), useMass_(false)
81{ }
82
83// Use the factory to build the preconditioner (this is where the work goes)
84LinearOp SIMPLEPreconditionerFactory
85 ::buildPreconditionerOperator(BlockedLinearOp & blockOp,
86 BlockPreconditionerState & state) const
87{
88 Teko_DEBUG_SCOPE("SIMPLEPreconditionerFactory::buildPreconditionerOperator",10);
89 Teko_DEBUG_EXPR(Teuchos::Time timer(""));
90
91 int rows = blockRowCount(blockOp);
92 int cols = blockColCount(blockOp);
93
94 TEUCHOS_ASSERT(rows==2); // sanity checks
95 TEUCHOS_ASSERT(cols==2);
96
97 bool buildExplicitSchurComplement = true;
98
99 // extract subblocks
100 const LinearOp F = getBlock(0,0,blockOp);
101 const LinearOp Bt = getBlock(0,1,blockOp);
102 const LinearOp B = getBlock(1,0,blockOp);
103 const LinearOp C = getBlock(1,1,blockOp);
104
105 LinearOp matF = F;
106 if(useMass_) {
107 TEUCHOS_ASSERT(massMatrix_!=Teuchos::null);
108 matF = massMatrix_;
109 }
110
111 // get approximation of inv(F) name H
112 std::string fApproxStr = "<error>";
113 LinearOp H;
114 if(fInverseType_==NotDiag) {
115 H = buildInverse(*customHFactory_,matF);
116 fApproxStr = customHFactory_->toString();
117
118 // since H is now implicit, we must build an implicit Schur complement
119 buildExplicitSchurComplement = false;
120 }
121 else if(fInverseType_==BlkDiag) {
122#ifdef TEKO_HAVE_EPETRA
123 // Block diagonal approximation for H
126 Hfact.initializeFromParameterList(BlkDiagList_);
127 H = Hfact.buildPreconditionerOperator(matF,Hstate);
128
129/*
130 // Get a FECrsMarix out of the BDP
131 RCP<Epetra_FECrsMatrix> Hcrs=rcp(Hstate.BDP_->CreateFECrsMatrix());
132 H=Thyra::epetraLinearOp(Hcrs);
133*/
134
135 buildExplicitSchurComplement = true; // NTS: Do I need this?
136 // Answer - no, but it is documenting whats going on here.
137#else
138 throw std::logic_error("SIMPLEPreconditionerFactory fInverseType_ == "
139 "BlkDiag but EPETRA is turned off!");
140#endif
141
142 }
143 else {
144 // get generic diagonal
145 H = getInvDiagonalOp(matF,fInverseType_);
146 fApproxStr = getDiagonalName(fInverseType_);
147 }
148
149 // adjust H for time scaling if it is a mass matrix
150 if(useMass_) {
151 RCP<const Teuchos::ParameterList> pl = state.getParameterList();
152
153 if(pl->isParameter("stepsize")) {
154 // get the step size
155 double stepsize = pl->get<double>("stepsize");
156
157 // scale by stepsize only if it is larger than 0
158 if(stepsize>0.0)
159 H = scale(stepsize,H);
160 }
161 }
162
163 // build approximate Schur complement: hatS = -C + B*H*Bt
164 LinearOp HBt, hatS;
165
166 if(buildExplicitSchurComplement) {
167 ModifiableLinearOp & mHBt = state.getModifiableOp("HBt");
168 ModifiableLinearOp & mhatS = state.getModifiableOp("hatS");
169 ModifiableLinearOp & BHBt = state.getModifiableOp("BHBt");
170
171 // build H*Bt
172 mHBt = explicitMultiply(H,Bt,mHBt);
173 HBt = mHBt;
174
175 // build B*H*Bt
176 BHBt = explicitMultiply(B,HBt,BHBt);
177
178 // build C-B*H*Bt
179 mhatS = explicitAdd(C,scale(-1.0,BHBt),mhatS);
180 hatS = mhatS;
181 }
182 else {
183 // build an implicit Schur complement
184 HBt = multiply(H,Bt);
185
186 hatS = add(C,scale(-1.0,multiply(B,HBt)));
187 }
188
189 // build the inverse for F
190 ModifiableLinearOp & invF = state.getModifiableOp("invF");
191 if(invF==Teuchos::null)
192 invF = buildInverse(*invVelFactory_,F);
193 else
194 rebuildInverse(*invVelFactory_,F,invF);
195
196 // build the approximate Schur complement
197 ModifiableLinearOp & invS = state.getModifiableOp("invS");
198 if(invS==Teuchos::null)
199 invS = buildInverse(*invPrsFactory_,hatS);
200 else
201 rebuildInverse(*invPrsFactory_,hatS,invS);
202
203 std::vector<LinearOp> invDiag(2); // vector storing inverses
204
205 // build lower triangular inverse matrix
206 BlockedLinearOp L = zeroBlockedOp(blockOp);
207 setBlock(1,0,L,B);
208 endBlockFill(L);
209
210 invDiag[0] = invF;
211 invDiag[1] = invS;
212 LinearOp invL = createBlockLowerTriInverseOp(L,invDiag);
213
214 // build upper triangular matrix
215 BlockedLinearOp U = zeroBlockedOp(blockOp);
216 setBlock(0,1,U,scale(1.0/alpha_,HBt));
217 endBlockFill(U);
218
219 invDiag[0] = identity(rangeSpace(invF));
220 invDiag[1] = scale(alpha_,identity(rangeSpace(invS)));
221 LinearOp invU = createBlockUpperTriInverseOp(U,invDiag);
222
223 // return implicit product operator
224 return multiply(invU,invL,"SIMPLE_"+fApproxStr);
225}
226
228void SIMPLEPreconditionerFactory::initializeFromParameterList(const Teuchos::ParameterList & pl)
229{
230 RCP<const InverseLibrary> invLib = getInverseLibrary();
231
232 // default conditions
233 useMass_ = false;
234 customHFactory_ = Teuchos::null;
235 fInverseType_ = Diagonal;
236
237 // get string specifying inverse
238 std::string invStr="", invVStr="", invPStr="";
239 alpha_ = 1.0;
240
241 // "parse" the parameter list
242 if(pl.isParameter("Inverse Type"))
243 invStr = pl.get<std::string>("Inverse Type");
244 if(pl.isParameter("Inverse Velocity Type"))
245 invVStr = pl.get<std::string>("Inverse Velocity Type");
246 if(pl.isParameter("Inverse Pressure Type"))
247 invPStr = pl.get<std::string>("Inverse Pressure Type");
248 if(pl.isParameter("Alpha"))
249 alpha_ = pl.get<double>("Alpha");
250 if(pl.isParameter("Explicit Velocity Inverse Type")) {
251 std::string fInverseStr = pl.get<std::string>("Explicit Velocity Inverse Type");
252
253 // build inverse types
254 fInverseType_ = getDiagonalType(fInverseStr);
255 if(fInverseType_==NotDiag)
256 customHFactory_ = invLib->getInverseFactory(fInverseStr);
257
258 // Grab the sublist if we're using the block diagonal
259 if(fInverseType_==BlkDiag)
260 BlkDiagList_=pl.sublist("H options");
261 }
262 if(pl.isParameter("Use Mass Scaling"))
263 useMass_ = pl.get<bool>("Use Mass Scaling");
264
265 Teko_DEBUG_MSG_BEGIN(5)
266 DEBUG_STREAM << "SIMPLE Parameters: " << std::endl;
267 DEBUG_STREAM << " inv type = \"" << invStr << "\"" << std::endl;
268 DEBUG_STREAM << " inv v type = \"" << invVStr << "\"" << std::endl;
269 DEBUG_STREAM << " inv p type = \"" << invPStr << "\"" << std::endl;
270 DEBUG_STREAM << " alpha = " << alpha_ << std::endl;
271 DEBUG_STREAM << " use mass = " << useMass_ << std::endl;
272 DEBUG_STREAM << " vel scaling = " << getDiagonalName(fInverseType_) << std::endl;
273 DEBUG_STREAM << "SIMPLE Parameter list: " << std::endl;
274 pl.print(DEBUG_STREAM);
275 Teko_DEBUG_MSG_END()
276
277 // set defaults as needed
278 if(invStr=="") invStr = "Amesos";
279 if(invVStr=="") invVStr = invStr;
280 if(invPStr=="") invPStr = invStr;
281
282 // two inverse factory objects
283 RCP<InverseFactory> invVFact, invPFact;
284
285 // build velocity inverse factory
286 invVFact = invLib->getInverseFactory(invVStr);
287 invPFact = invVFact; // by default these are the same
288 if(invVStr!=invPStr) // if different, build pressure inverse factory
289 invPFact = invLib->getInverseFactory(invPStr);
290
291 // based on parameter type build a strategy
292 invVelFactory_ = invVFact;
293 invPrsFactory_ = invPFact;
294
295 if(useMass_) {
296 Teuchos::RCP<Teko::RequestHandler> rh = getRequestHandler();
297 rh->preRequest<Teko::LinearOp>(Teko::RequestMesg("Velocity Mass Matrix"));
298 Teko::LinearOp mass
299 = rh->request<Teko::LinearOp>(Teko::RequestMesg("Velocity Mass Matrix"));
300 setMassMatrix(mass);
301 }
302}
303
305Teuchos::RCP<Teuchos::ParameterList> SIMPLEPreconditionerFactory::getRequestedParameters() const
306{
307 Teuchos::RCP<Teuchos::ParameterList> result;
308 Teuchos::RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList());
309
310 // grab parameters from F solver
311 RCP<Teuchos::ParameterList> vList = invVelFactory_->getRequestedParameters();
312 if(vList!=Teuchos::null) {
313 Teuchos::ParameterList::ConstIterator itr;
314 for(itr=vList->begin();itr!=vList->end();++itr)
315 pl->setEntry(itr->first,itr->second);
316 result = pl;
317 }
318
319 // grab parameters from S solver
320 RCP<Teuchos::ParameterList> pList = invPrsFactory_->getRequestedParameters();
321 if(pList!=Teuchos::null) {
322 Teuchos::ParameterList::ConstIterator itr;
323 for(itr=pList->begin();itr!=pList->end();++itr)
324 pl->setEntry(itr->first,itr->second);
325 result = pl;
326 }
327
328 // grab parameters from S solver
329 if(customHFactory_!=Teuchos::null) {
330 RCP<Teuchos::ParameterList> hList = customHFactory_->getRequestedParameters();
331 if(hList!=Teuchos::null) {
332 Teuchos::ParameterList::ConstIterator itr;
333 for(itr=hList->begin();itr!=hList->end();++itr)
334 pl->setEntry(itr->first,itr->second);
335 result = pl;
336 }
337 }
338
339 return result;
340}
341
343bool SIMPLEPreconditionerFactory::updateRequestedParameters(const Teuchos::ParameterList & pl)
344{
345 Teko_DEBUG_SCOPE("InvLSCStrategy::updateRequestedParameters",10);
346 bool result = true;
347
348 // update requested parameters in solvers
349 result &= invVelFactory_->updateRequestedParameters(pl);
350 result &= invPrsFactory_->updateRequestedParameters(pl);
351 if(customHFactory_!=Teuchos::null)
352 result &= customHFactory_->updateRequestedParameters(pl);
353
354 return result;
355}
356
357} // end namespace NS
358} // end namespace Teko
InverseLinearOp buildInverse(const InverseFactory &factory, const LinearOp &A)
Build an inverse operator using a factory and a linear operator.
void rebuildInverse(const InverseFactory &factory, const LinearOp &A, InverseLinearOp &invA)
void scale(const double alpha, MultiVector &x)
Scale a multivector by a constant.
@ BlkDiag
Specifies that a block diagonal approximation is to be used.
@ NotDiag
For user convenience, if Teko recieves this value, exceptions will be thrown.
@ Diagonal
Specifies that just the diagonal is used.
int blockRowCount(const BlockedLinearOp &blo)
Get the row count in a block linear operator.
MultiVector getBlock(int i, const BlockedMultiVector &bmv)
Get the ith block from a BlockedMultiVector object.
void setBlock(int i, int j, BlockedLinearOp &blo, const LinearOp &lo)
Set the i,j block in a BlockedLinearOp object.
VectorSpace rangeSpace(const LinearOp &lo)
Replace nonzeros with a scalar value, used to zero out an operator.
int blockColCount(const BlockedLinearOp &blo)
Get the column count in a block linear operator.
void endBlockFill(BlockedLinearOp &blo)
Notify the blocked operator that the fill stage is completed.
An implementation of a state object for block preconditioners.
Preconditioner factory for building explcit inverse of diagonal operators. This includes block operat...
virtual void initializeFromParameterList(const Teuchos::ParameterList &pl)
Initialize from a parameter list.
LinearOp buildPreconditionerOperator(LinearOp &lo, PreconditionerState &state) const
virtual bool updateRequestedParameters(const Teuchos::ParameterList &pl)
For assisting in construction of the preconditioner.
virtual void initializeFromParameterList(const Teuchos::ParameterList &pl)
Initialize from a parameter list.
virtual void setMassMatrix(Teko::LinearOp &mass)
Set the mass matrix for this factory.
virtual Teuchos::RCP< Teuchos::ParameterList > getRequestedParameters() const
For assisting in construction of the preconditioner.
Teuchos::RCP< const InverseLibrary > getInverseLibrary() const
Get the inverse library used by this preconditioner factory.
Teuchos::RCP< RequestHandler > getRequestHandler() const
Get the request handler with pointers to the appropriate callbacks.
virtual Teko::ModifiableLinearOp & getModifiableOp(const std::string &name)
Add a named operator to the state object.