Teko Version of the Day
Loading...
Searching...
No Matches
Teko_TimingsSIMPLEPreconditionerFactory.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_TimingsSIMPLEPreconditionerFactory.hpp"
48
49#include "Teko_Utilities.hpp"
50#include "Teko_InverseFactory.hpp"
51#include "Teko_BlockLowerTriInverseOp.hpp"
52#include "Teko_BlockUpperTriInverseOp.hpp"
53#ifdef TEKO_HAVE_EPETRA
54#include "Teko_DiagonalPreconditionerFactory.hpp"
55#endif
56
57#include "Teuchos_Time.hpp"
58
59using Teuchos::RCP;
60
61namespace Teko {
62namespace NS {
63
64// Constructor definition
65TimingsSIMPLEPreconditionerFactory
66 ::TimingsSIMPLEPreconditionerFactory(const RCP<InverseFactory> & inverse,
67 double alpha)
68 : invVelFactory_(inverse), invPrsFactory_(inverse), alpha_(alpha), fInverseType_(Diagonal), useMass_(false)
69 , constrTotal_("SIMPLE Constr: Total"), subTotal_("SIMPLE Constr: Subs"), constrCount_(0)
70{ }
71
72TimingsSIMPLEPreconditionerFactory
73 ::TimingsSIMPLEPreconditionerFactory(const RCP<InverseFactory> & invVFact,
74 const RCP<InverseFactory> & invPFact,
75 double alpha)
76 : invVelFactory_(invVFact), invPrsFactory_(invPFact), alpha_(alpha), fInverseType_(Diagonal), useMass_(false)
77 , constrTotal_("SIMPLE Constr: Total"), subTotal_("SIMPLE Constr: Subs"), constrCount_(0)
78{ }
79
80TimingsSIMPLEPreconditionerFactory::TimingsSIMPLEPreconditionerFactory()
81 : alpha_(1.0), fInverseType_(Diagonal), useMass_(false)
82 , constrTotal_("SIMPLE Constr: Total"), subTotal_("SIMPLE Constr: Subs"), constrCount_(0)
83{ }
84
86{
87 if(constrTotal_.totalElapsedTime()>0.0) {
88 Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
89 out.setOutputToRootOnly(0);
90
91 out << "===========================================================================" << std::endl;
92 out << std::endl;
93 out << "SIMPLE Construction Count = " << constrCount_ << std::endl;
94 out << "SIMPLE Construction Total = " << constrTotal_.totalElapsedTime() << std::endl;
95 out << "SIMPLE Sub Components Total = " << subTotal_.totalElapsedTime() << std::endl;
96 out << std::endl;
97 out << "===========================================================================" << std::endl;
98 }
99}
100
101// Use the factory to build the preconditioner (this is where the work goes)
102LinearOp TimingsSIMPLEPreconditionerFactory
103 ::buildPreconditionerOperator(BlockedLinearOp & blockOp,
104 BlockPreconditionerState & state) const
105{
106 constrTotal_.start();
107
108 Teko_DEBUG_SCOPE("TimingsSIMPLEPreconditionerFactory::buildPreconditionerOperator",10);
109 Teko_DEBUG_EXPR(Teuchos::Time timer(""));
110
111 int rows = blockRowCount(blockOp);
112 int cols = blockColCount(blockOp);
113
114 TEUCHOS_ASSERT(rows==2); // sanity checks
115 TEUCHOS_ASSERT(cols==2);
116
117 bool buildExplicitSchurComplement = true;
118
119 // extract subblocks
120 const LinearOp F = getBlock(0,0,blockOp);
121 const LinearOp Bt = getBlock(0,1,blockOp);
122 const LinearOp B = getBlock(1,0,blockOp);
123 const LinearOp C = getBlock(1,1,blockOp);
124
125 LinearOp matF = F;
126 if(useMass_) {
127 TEUCHOS_ASSERT(massMatrix_!=Teuchos::null);
128 matF = massMatrix_;
129 }
130
131 // get approximation of inv(F) name H
132 std::string fApproxStr = "<error>";
133 LinearOp H;
134 if(fInverseType_==NotDiag) {
135 H = buildInverse(*customHFactory_,matF);
136 fApproxStr = customHFactory_->toString();
137
138 // since H is now implicit, we must build an implicit Schur complement
139 buildExplicitSchurComplement = false;
140 }
141 else if(fInverseType_==BlkDiag) {
142#ifdef TEKO_HAVE_EPETRA
143 // Block diagonal approximation for H
146 Hfact.initializeFromParameterList(BlkDiagList_);
147 H = Hfact.buildPreconditionerOperator(matF,Hstate);
148
149 buildExplicitSchurComplement = true; // NTS: Do I need this?
150 // Answer - no, but it is documenting whats going on here.
151#else
152 throw std::logic_error("TimingsSIMPLEPreconditionerFactory fInverseType_ == "
153 "BlkDiag but EPETRA is turned off!");
154#endif
155 }
156 else {
157 // get generic diagonal
158 subTotal_.start();
159 H = getInvDiagonalOp(matF,fInverseType_);
160 subTotal_.stop();
161 fApproxStr = getDiagonalName(fInverseType_);
162 }
163
164 // adjust H for time scaling if it is a mass matrix
165 if(useMass_) {
166 RCP<const Teuchos::ParameterList> pl = state.getParameterList();
167
168 if(pl->isParameter("stepsize")) {
169 // get the step size
170 double stepsize = pl->get<double>("stepsize");
171
172 // scale by stepsize only if it is larger than 0
173 if(stepsize>0.0)
174 H = scale(stepsize,H);
175 }
176 }
177
178 // build approximate Schur complement: hatS = -C + B*H*Bt
179 LinearOp HBt, hatS;
180
181 if(buildExplicitSchurComplement) {
182 ModifiableLinearOp & mHBt = state.getModifiableOp("HBt");
183 ModifiableLinearOp & mhatS = state.getModifiableOp("hatS");
184 ModifiableLinearOp & BHBt = state.getModifiableOp("BHBt");
185
186 // build H*Bt
187 subTotal_.start();
188 mHBt = explicitMultiply(H,Bt,mHBt);
189 subTotal_.stop();
190 HBt = mHBt;
191
192 // build B*H*Bt
193 subTotal_.start();
194 BHBt = explicitMultiply(B,HBt,BHBt);
195 subTotal_.stop();
196
197 // build C-B*H*Bt
198 subTotal_.start();
199 mhatS = explicitAdd(C,scale(-1.0,BHBt),mhatS);
200 subTotal_.stop();
201 hatS = mhatS;
202 }
203 else {
204 // build an implicit Schur complement
205 HBt = multiply(H,Bt);
206
207 hatS = add(C,scale(-1.0,multiply(B,HBt)));
208 }
209
210 // time the application of HBt
211 if(timed_HBt_==Teuchos::null) {
212 timed_HBt_ = Teuchos::rcp(new DiagnosticLinearOp(getOutputStream(),HBt,"HBt"));
213 }
214 else {
215 timed_HBt_->setLinearOp(HBt);
216 }
217
218 // time the application of B
219 if(timed_B_==Teuchos::null) {
220 timed_B_ = Teuchos::rcp(new DiagnosticLinearOp(getOutputStream(),B,"B"));
221 }
222 else {
223 timed_B_->setLinearOp(B);
224 }
225
226 // build the inverse for F
227 ModifiableLinearOp & invF = state.getModifiableOp("invF");
228 subTotal_.start();
229 if(invF==Teuchos::null) {
230 invF = buildInverse(*invVelFactory_,F);
231
232 timed_invF_ = Teuchos::rcp(new DiagnosticLinearOp(getOutputStream(),invF,"invF"));
233 }
234 else {
235 rebuildInverse(*invVelFactory_,F,invF);
236
237 timed_invF_->setLinearOp(invF);
238 }
239 subTotal_.stop();
240
241 // build the approximate Schur complement
242 ModifiableLinearOp & invS = state.getModifiableOp("invS");
243 subTotal_.start();
244 if(invS==Teuchos::null) {
245 invS = buildInverse(*invPrsFactory_,hatS);
246
247 timed_invS_ = Teuchos::rcp(new DiagnosticLinearOp(getOutputStream(),invS,"invS"));
248 }
249 else {
250 rebuildInverse(*invPrsFactory_,hatS,invS);
251
252 timed_invS_->setLinearOp(invS);
253 }
254 subTotal_.stop();
255
256 std::vector<LinearOp> invDiag(2); // vector storing inverses
257
258 // build lower triangular inverse matrix
259 BlockedLinearOp L = zeroBlockedOp(blockOp);
260 setBlock(1,0,L,timed_B_);
261 endBlockFill(L);
262
263 invDiag[0] = timed_invF_;
264 invDiag[1] = timed_invS_;
265 LinearOp invL = createBlockLowerTriInverseOp(L,invDiag);
266
267 // build upper triangular matrix
268 BlockedLinearOp U = zeroBlockedOp(blockOp);
269 setBlock(0,1,U,scale<double>(1.0/alpha_,timed_HBt_.getConst()));
270 endBlockFill(U);
271
272 invDiag[0] = identity(rangeSpace(invF));
273 invDiag[1] = scale(alpha_,identity(rangeSpace(invS)));
274 LinearOp invU = createBlockUpperTriInverseOp(U,invDiag);
275
276 // return implicit product operator
277 Teko::LinearOp iU_t_iL = multiply(invU,invL,"SIMPLE_"+fApproxStr);
278
279 // time the application of iU_t_iL
280 if(timed_iU_t_iL_==Teuchos::null)
281 timed_iU_t_iL_ = Teuchos::rcp(new DiagnosticLinearOp(getOutputStream(),iU_t_iL,"iU_t_iL"));
282 else
283 timed_iU_t_iL_->setLinearOp(iU_t_iL);
284
285 constrCount_++;
286
287 constrTotal_.stop();
288
289 return timed_iU_t_iL_;
290}
291
294{
295 RCP<const InverseLibrary> invLib = getInverseLibrary();
296
297 // default conditions
298 useMass_ = false;
299 customHFactory_ = Teuchos::null;
300 fInverseType_ = Diagonal;
301
302 // get string specifying inverse
303 std::string invStr="", invVStr="", invPStr="";
304 alpha_ = 1.0;
305
306 // "parse" the parameter list
307 if(pl.isParameter("Inverse Type"))
308 invStr = pl.get<std::string>("Inverse Type");
309 if(pl.isParameter("Inverse Velocity Type"))
310 invVStr = pl.get<std::string>("Inverse Velocity Type");
311 if(pl.isParameter("Inverse Pressure Type"))
312 invPStr = pl.get<std::string>("Inverse Pressure Type");
313 if(pl.isParameter("Alpha"))
314 alpha_ = pl.get<double>("Alpha");
315 if(pl.isParameter("Explicit Velocity Inverse Type")) {
316 std::string fInverseStr = pl.get<std::string>("Explicit Velocity Inverse Type");
317
318 // build inverse types
319 fInverseType_ = getDiagonalType(fInverseStr);
320 if(fInverseType_==NotDiag)
321 customHFactory_ = invLib->getInverseFactory(fInverseStr);
322
323 // Grab the sublist if we're using the block diagonal
324 if(fInverseType_==BlkDiag)
325 BlkDiagList_=pl.sublist("H options");
326 }
327 if(pl.isParameter("Use Mass Scaling"))
328 useMass_ = pl.get<bool>("Use Mass Scaling");
329
330 Teko_DEBUG_MSG_BEGIN(5)
331 DEBUG_STREAM << "SIMPLE Parameters: " << std::endl;
332 DEBUG_STREAM << " inv type = \"" << invStr << "\"" << std::endl;
333 DEBUG_STREAM << " inv v type = \"" << invVStr << "\"" << std::endl;
334 DEBUG_STREAM << " inv p type = \"" << invPStr << "\"" << std::endl;
335 DEBUG_STREAM << " alpha = " << alpha_ << std::endl;
336 DEBUG_STREAM << " use mass = " << useMass_ << std::endl;
337 DEBUG_STREAM << " vel scaling = " << getDiagonalName(fInverseType_) << std::endl;
338 DEBUG_STREAM << "SIMPLE Parameter list: " << std::endl;
339 pl.print(DEBUG_STREAM);
340 Teko_DEBUG_MSG_END()
341
342 // set defaults as needed
343 if(invStr=="") invStr = "Amesos";
344 if(invVStr=="") invVStr = invStr;
345 if(invPStr=="") invPStr = invStr;
346
347 // two inverse factory objects
348 RCP<InverseFactory> invVFact, invPFact;
349
350 // build velocity inverse factory
351 invVFact = invLib->getInverseFactory(invVStr);
352 invPFact = invVFact; // by default these are the same
353 if(invVStr!=invPStr) // if different, build pressure inverse factory
354 invPFact = invLib->getInverseFactory(invPStr);
355
356 // based on parameter type build a strategy
357 invVelFactory_ = invVFact;
358 invPrsFactory_ = invPFact;
359
360 if(useMass_) {
361 Teuchos::RCP<Teko::RequestHandler> rh = getRequestHandler();
362 rh->preRequest<Teko::LinearOp>(Teko::RequestMesg("Velocity Mass Matrix"));
363 Teko::LinearOp mass
364 = rh->request<Teko::LinearOp>(Teko::RequestMesg("Velocity Mass Matrix"));
365 setMassMatrix(mass);
366 }
367}
368
370Teuchos::RCP<Teuchos::ParameterList> TimingsSIMPLEPreconditionerFactory::getRequestedParameters() const
371{
372 Teuchos::RCP<Teuchos::ParameterList> result;
373 Teuchos::RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList());
374
375 // grab parameters from F solver
376 RCP<Teuchos::ParameterList> vList = invVelFactory_->getRequestedParameters();
377 if(vList!=Teuchos::null) {
378 Teuchos::ParameterList::ConstIterator itr;
379 for(itr=vList->begin();itr!=vList->end();++itr)
380 pl->setEntry(itr->first,itr->second);
381 result = pl;
382 }
383
384 // grab parameters from S solver
385 RCP<Teuchos::ParameterList> pList = invPrsFactory_->getRequestedParameters();
386 if(pList!=Teuchos::null) {
387 Teuchos::ParameterList::ConstIterator itr;
388 for(itr=pList->begin();itr!=pList->end();++itr)
389 pl->setEntry(itr->first,itr->second);
390 result = pl;
391 }
392
393 // grab parameters from S solver
394 if(customHFactory_!=Teuchos::null) {
395 RCP<Teuchos::ParameterList> hList = customHFactory_->getRequestedParameters();
396 if(hList!=Teuchos::null) {
397 Teuchos::ParameterList::ConstIterator itr;
398 for(itr=hList->begin();itr!=hList->end();++itr)
399 pl->setEntry(itr->first,itr->second);
400 result = pl;
401 }
402 }
403
404 return result;
405}
406
409{
410 Teko_DEBUG_SCOPE("InvLSCStrategy::updateRequestedParameters",10);
411 bool result = true;
412
413 // update requested parameters in solvers
414 result &= invVelFactory_->updateRequestedParameters(pl);
415 result &= invPrsFactory_->updateRequestedParameters(pl);
416 if(customHFactory_!=Teuchos::null)
417 result &= customHFactory_->updateRequestedParameters(pl);
418
419 return result;
420}
421
422} // end namespace NS
423} // 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.
This linear operator prints diagnostics about operator application and creation times....
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 ~TimingsSIMPLEPreconditionerFactory()
Destructor that outputs construction timings.
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.