Teko  Version of the Day
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 #include "Teko_DiagonalPreconditionerFactory.hpp"
54 
55 #include "Teuchos_Time.hpp"
56 #include "Epetra_FECrsGraph.h"
57 #include "Epetra_FECrsMatrix.h"
58 #include "EpetraExt_PointToBlockDiagPermute.h"
59 #include "EpetraExt_MatrixMatrix.h"
60 #include "Thyra_EpetraOperatorWrapper.hpp"
61 #include "Thyra_EpetraLinearOp.hpp"
62 
63 
64 using Teuchos::RCP;
65 
66 namespace Teko {
67 namespace NS {
68 
69 // Constructor definition
71  ::TimingsSIMPLEPreconditionerFactory(const RCP<InverseFactory> & inverse,
72  double alpha)
73  : invVelFactory_(inverse), invPrsFactory_(inverse), alpha_(alpha), fInverseType_(Diagonal), useMass_(false)
74  , constrTotal_("SIMPLE Constr: Total"), subTotal_("SIMPLE Constr: Subs"), constrCount_(0)
75 { }
76 
78  ::TimingsSIMPLEPreconditionerFactory(const RCP<InverseFactory> & invVFact,
79  const RCP<InverseFactory> & invPFact,
80  double alpha)
81  : invVelFactory_(invVFact), invPrsFactory_(invPFact), alpha_(alpha), fInverseType_(Diagonal), useMass_(false)
82  , constrTotal_("SIMPLE Constr: Total"), subTotal_("SIMPLE Constr: Subs"), constrCount_(0)
83 { }
84 
86  : alpha_(1.0), fInverseType_(Diagonal), useMass_(false)
87  , constrTotal_("SIMPLE Constr: Total"), subTotal_("SIMPLE Constr: Subs"), constrCount_(0)
88 { }
89 
91 {
92  if(constrTotal_.totalElapsedTime()>0.0) {
93  Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
94  out.setOutputToRootOnly(0);
95 
96  out << "===========================================================================" << std::endl;
97  out << std::endl;
98  out << "SIMPLE Construction Count = " << constrCount_ << std::endl;
99  out << "SIMPLE Construction Total = " << constrTotal_.totalElapsedTime() << std::endl;
100  out << "SIMPLE Sub Components Total = " << subTotal_.totalElapsedTime() << std::endl;
101  out << std::endl;
102  out << "===========================================================================" << std::endl;
103  }
104 }
105 
106 // Use the factory to build the preconditioner (this is where the work goes)
108  ::buildPreconditionerOperator(BlockedLinearOp & blockOp,
109  BlockPreconditionerState & state) const
110 {
111  constrTotal_.start();
112 
113  Teko_DEBUG_SCOPE("TimingsSIMPLEPreconditionerFactory::buildPreconditionerOperator",10);
114  Teko_DEBUG_EXPR(Teuchos::Time timer(""));
115 
116  int rows = blockRowCount(blockOp);
117  int cols = blockColCount(blockOp);
118 
119  TEUCHOS_ASSERT(rows==2); // sanity checks
120  TEUCHOS_ASSERT(cols==2);
121 
122  bool buildExplicitSchurComplement = true;
123 
124  // extract subblocks
125  const LinearOp F = getBlock(0,0,blockOp);
126  const LinearOp Bt = getBlock(0,1,blockOp);
127  const LinearOp B = getBlock(1,0,blockOp);
128  const LinearOp C = getBlock(1,1,blockOp);
129 
130  LinearOp matF = F;
131  if(useMass_) {
132  TEUCHOS_ASSERT(massMatrix_!=Teuchos::null);
133  matF = massMatrix_;
134  }
135 
136  // get approximation of inv(F) name H
137  std::string fApproxStr = "<error>";
138  LinearOp H;
139  if(fInverseType_==NotDiag) {
140  H = buildInverse(*customHFactory_,matF);
141  fApproxStr = customHFactory_->toString();
142 
143  // since H is now implicit, we must build an implicit Schur complement
144  buildExplicitSchurComplement = false;
145  }
146  else if(fInverseType_==BlkDiag) {
147  // Block diagonal approximation for H
149  DiagonalPrecondState Hstate;
150  Hfact.initializeFromParameterList(BlkDiagList_);
151  H = Hfact.buildPreconditionerOperator(matF,Hstate);
152 
153  buildExplicitSchurComplement = true; // NTS: Do I need this?
154  // Answer - no, but it is documenting whats going on here.
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 
370 Teuchos::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
virtual bool updateRequestedParameters(const Teuchos::ParameterList &pl)
For assisting in construction of the preconditioner.
virtual ~TimingsSIMPLEPreconditionerFactory()
Destructor that outputs construction timings.
Specifies that a block diagonal approximation is to be used.
Teuchos::RCP< RequestHandler > getRequestHandler() const
Get the request handler with pointers to the appropriate callbacks.
Specifies that just the diagonal is used.
int blockRowCount(const BlockedLinearOp &blo)
Get the row count in a block linear operator.
InverseLinearOp buildInverse(const InverseFactory &factory, const LinearOp &A)
Build an inverse operator using a factory and a linear operator.
virtual void setMassMatrix(Teko::LinearOp &mass)
Set the mass matrix for this factory.
Teuchos::RCP< const InverseLibrary > getInverseLibrary() const
Get the inverse library used by this preconditioner factory.
void scale(const double alpha, MultiVector &x)
Scale a multivector by a constant.
LinearOp buildPreconditionerOperator(BlockedLinearOp &blo, BlockPreconditionerState &state) const
Function that is called to build the preconditioner for the linear operator that is passed in...
LinearOp buildPreconditionerOperator(LinearOp &lo, PreconditionerState &state) const
void endBlockFill(BlockedLinearOp &blo)
Notify the blocked operator that the fill stage is completed.
virtual void initializeFromParameterList(const Teuchos::ParameterList &pl)
Initialize from a parameter list.
Preconditioner factory for building explcit inverse of diagonal operators. This includes block operat...
void rebuildInverse(const InverseFactory &factory, const LinearOp &A, InverseLinearOp &invA)
An implementation of a state object for block preconditioners.
VectorSpace rangeSpace(const LinearOp &lo)
Get the range space of a linear operator.
virtual Teuchos::RCP< Teuchos::ParameterList > getRequestedParameters() const
For assisting in construction of the preconditioner.
This linear operator prints diagnostics about operator application and creation times. It is useful for debugging problems and determining bottle necks.
MultiVector getBlock(int i, const BlockedMultiVector &bmv)
Get the ith block from a BlockedMultiVector object.
For user convenience, if Teko recieves this value, exceptions will be thrown.
virtual void initializeFromParameterList(const Teuchos::ParameterList &pl)
Initialize from a parameter list.
int blockColCount(const BlockedLinearOp &blo)
Get the column count in a block linear operator.
void setBlock(int i, int j, BlockedLinearOp &blo, const LinearOp &lo)
Set the i,j block in a BlockedLinearOp object.
virtual Teko::ModifiableLinearOp & getModifiableOp(const std::string &name)
Add a named operator to the state object.