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