Teko  Version of the Day
 All Classes Files Functions Variables Pages
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 
60 using Teuchos::RCP;
61 
62 namespace Teko {
63 namespace NS {
64 
65 // Constructor definition
66 SIMPLEPreconditionerFactory ::SIMPLEPreconditionerFactory(const RCP<InverseFactory>& inverse,
67  double alpha)
68  : invVelFactory_(inverse),
69  invPrsFactory_(inverse),
70  alpha_(alpha),
71  fInverseType_(Diagonal),
72  useMass_(false) {}
73 
74 SIMPLEPreconditionerFactory ::SIMPLEPreconditionerFactory(const RCP<InverseFactory>& invVFact,
75  const RCP<InverseFactory>& invPFact,
76  double alpha)
77  : invVelFactory_(invVFact),
78  invPrsFactory_(invPFact),
79  alpha_(alpha),
80  fInverseType_(Diagonal),
81  useMass_(false) {}
82 
83 SIMPLEPreconditionerFactory::SIMPLEPreconditionerFactory()
84  : alpha_(1.0), fInverseType_(Diagonal), useMass_(false) {}
85 
86 // Use the factory to build the preconditioner (this is where the work goes)
88  BlockedLinearOp& blockOp, BlockPreconditionerState& state) const {
89  Teko_DEBUG_SCOPE("SIMPLEPreconditionerFactory::buildPreconditionerOperator", 10);
90  Teko_DEBUG_EXPR(Teuchos::Time timer(""));
91 
92  int rows = blockRowCount(blockOp);
93  int cols = blockColCount(blockOp);
94 
95  TEUCHOS_ASSERT(rows == 2); // sanity checks
96  TEUCHOS_ASSERT(cols == 2);
97 
98  bool buildExplicitSchurComplement = true;
99 
100  // extract subblocks
101  const LinearOp F = getBlock(0, 0, blockOp);
102  const LinearOp Bt = getBlock(0, 1, blockOp);
103  const LinearOp B = getBlock(1, 0, blockOp);
104  const LinearOp C = getBlock(1, 1, blockOp);
105 
106  LinearOp matF = F;
107  if (useMass_) {
108  TEUCHOS_ASSERT(massMatrix_ != Teuchos::null);
109  matF = massMatrix_;
110  }
111 
112  // get approximation of inv(F) name H
113  std::string fApproxStr = "<error>";
114  LinearOp H;
115  if (fInverseType_ == NotDiag) {
116  H = buildInverse(*customHFactory_, matF);
117  fApproxStr = customHFactory_->toString();
118 
119  // since H is now implicit, we must build an implicit Schur complement
120  buildExplicitSchurComplement = false;
121  } else if (fInverseType_ == BlkDiag) {
122 #ifdef TEKO_HAVE_EPETRA
123  // Block diagonal approximation for H
125  DiagonalPrecondState Hstate;
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(
139  "SIMPLEPreconditionerFactory fInverseType_ == "
140  "BlkDiag but EPETRA is turned off!");
141 #endif
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) H = scale(stepsize, H);
159  }
160  }
161 
162  // build approximate Schur complement: hatS = -C + B*H*Bt
163  LinearOp HBt, hatS;
164 
165  if (buildExplicitSchurComplement) {
166  ModifiableLinearOp& mHBt = state.getModifiableOp("HBt");
167  ModifiableLinearOp& mhatS = state.getModifiableOp("hatS");
168  ModifiableLinearOp& BHBt = state.getModifiableOp("BHBt");
169 
170  // build H*Bt
171  mHBt = explicitMultiply(H, Bt, mHBt);
172  HBt = mHBt;
173 
174  // build B*H*Bt
175  BHBt = explicitMultiply(B, HBt, BHBt);
176 
177  // build C-B*H*Bt
178  mhatS = explicitAdd(C, scale(-1.0, BHBt), mhatS);
179  hatS = mhatS;
180  } else {
181  // build an implicit Schur complement
182  HBt = multiply(H, Bt);
183 
184  hatS = add(C, scale(-1.0, multiply(B, HBt)));
185  }
186 
187  // build the inverse for F
188  ModifiableLinearOp& invF = state.getModifiableOp("invF");
189  if (invF == Teuchos::null)
190  invF = buildInverse(*invVelFactory_, F);
191  else
192  rebuildInverse(*invVelFactory_, F, invF);
193 
194  // build the approximate Schur complement
195  ModifiableLinearOp& invS = state.getModifiableOp("invS");
196  if (invS == Teuchos::null)
197  invS = buildInverse(*invPrsFactory_, hatS);
198  else
199  rebuildInverse(*invPrsFactory_, hatS, invS);
200 
201  std::vector<LinearOp> invDiag(2); // vector storing inverses
202 
203  // build lower triangular inverse matrix
204  BlockedLinearOp L = zeroBlockedOp(blockOp);
205  setBlock(1, 0, L, B);
206  endBlockFill(L);
207 
208  invDiag[0] = invF;
209  invDiag[1] = invS;
210  LinearOp invL = createBlockLowerTriInverseOp(L, invDiag);
211 
212  // build upper triangular matrix
213  BlockedLinearOp U = zeroBlockedOp(blockOp);
214  setBlock(0, 1, U, scale(1.0 / alpha_, HBt));
215  endBlockFill(U);
216 
217  invDiag[0] = identity(rangeSpace(invF));
218  invDiag[1] = scale(alpha_, identity(rangeSpace(invS)));
219  LinearOp invU = createBlockUpperTriInverseOp(U, invDiag);
220 
221  // return implicit product operator
222  return multiply(invU, invL, "SIMPLE_" + fApproxStr);
223 }
224 
226 void SIMPLEPreconditionerFactory::initializeFromParameterList(const Teuchos::ParameterList& pl) {
227  RCP<const InverseLibrary> invLib = getInverseLibrary();
228 
229  // default conditions
230  useMass_ = false;
231  customHFactory_ = Teuchos::null;
232  fInverseType_ = Diagonal;
233 
234  // get string specifying inverse
235  std::string invStr = "", invVStr = "", invPStr = "";
236  alpha_ = 1.0;
237 
238  // "parse" the parameter list
239  if (pl.isParameter("Inverse Type")) invStr = pl.get<std::string>("Inverse Type");
240  if (pl.isParameter("Inverse Velocity Type"))
241  invVStr = pl.get<std::string>("Inverse Velocity Type");
242  if (pl.isParameter("Inverse Pressure Type"))
243  invPStr = pl.get<std::string>("Inverse Pressure Type");
244  if (pl.isParameter("Alpha")) alpha_ = pl.get<double>("Alpha");
245  if (pl.isParameter("Explicit Velocity Inverse Type")) {
246  std::string fInverseStr = pl.get<std::string>("Explicit Velocity Inverse Type");
247 
248  // build inverse types
249  fInverseType_ = getDiagonalType(fInverseStr);
250  if (fInverseType_ == NotDiag) customHFactory_ = invLib->getInverseFactory(fInverseStr);
251 
252  // Grab the sublist if we're using the block diagonal
253  if (fInverseType_ == BlkDiag) BlkDiagList_ = pl.sublist("H options");
254  }
255  if (pl.isParameter("Use Mass Scaling")) useMass_ = pl.get<bool>("Use Mass Scaling");
256 
257  Teko_DEBUG_MSG_BEGIN(5) DEBUG_STREAM << "SIMPLE Parameters: " << std::endl;
258  DEBUG_STREAM << " inv type = \"" << invStr << "\"" << std::endl;
259  DEBUG_STREAM << " inv v type = \"" << invVStr << "\"" << std::endl;
260  DEBUG_STREAM << " inv p type = \"" << invPStr << "\"" << std::endl;
261  DEBUG_STREAM << " alpha = " << alpha_ << std::endl;
262  DEBUG_STREAM << " use mass = " << useMass_ << std::endl;
263  DEBUG_STREAM << " vel scaling = " << getDiagonalName(fInverseType_) << std::endl;
264  DEBUG_STREAM << "SIMPLE Parameter list: " << std::endl;
265  pl.print(DEBUG_STREAM);
266  Teko_DEBUG_MSG_END()
267 
268  // set defaults as needed
269  if (invStr == "") invStr = "Amesos";
270  if (invVStr == "") invVStr = invStr;
271  if (invPStr == "") invPStr = invStr;
272 
273  // two inverse factory objects
274  RCP<InverseFactory> invVFact, invPFact;
275 
276  // build velocity inverse factory
277  invVFact = invLib->getInverseFactory(invVStr);
278  invPFact = invVFact; // by default these are the same
279  if (invVStr != invPStr) // if different, build pressure inverse factory
280  invPFact = invLib->getInverseFactory(invPStr);
281 
282  // based on parameter type build a strategy
283  invVelFactory_ = invVFact;
284  invPrsFactory_ = invPFact;
285 
286  if (useMass_) {
287  Teuchos::RCP<Teko::RequestHandler> rh = getRequestHandler();
288  rh->preRequest<Teko::LinearOp>(Teko::RequestMesg("Velocity Mass Matrix"));
289  Teko::LinearOp mass = rh->request<Teko::LinearOp>(Teko::RequestMesg("Velocity Mass Matrix"));
290  setMassMatrix(mass);
291  }
292 }
293 
295 Teuchos::RCP<Teuchos::ParameterList> SIMPLEPreconditionerFactory::getRequestedParameters() const {
296  Teuchos::RCP<Teuchos::ParameterList> result;
297  Teuchos::RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList());
298 
299  // grab parameters from F solver
300  RCP<Teuchos::ParameterList> vList = invVelFactory_->getRequestedParameters();
301  if (vList != Teuchos::null) {
302  Teuchos::ParameterList::ConstIterator itr;
303  for (itr = vList->begin(); itr != vList->end(); ++itr) pl->setEntry(itr->first, itr->second);
304  result = pl;
305  }
306 
307  // grab parameters from S solver
308  RCP<Teuchos::ParameterList> pList = invPrsFactory_->getRequestedParameters();
309  if (pList != Teuchos::null) {
310  Teuchos::ParameterList::ConstIterator itr;
311  for (itr = pList->begin(); itr != pList->end(); ++itr) pl->setEntry(itr->first, itr->second);
312  result = pl;
313  }
314 
315  // grab parameters from S solver
316  if (customHFactory_ != Teuchos::null) {
317  RCP<Teuchos::ParameterList> hList = customHFactory_->getRequestedParameters();
318  if (hList != Teuchos::null) {
319  Teuchos::ParameterList::ConstIterator itr;
320  for (itr = hList->begin(); itr != hList->end(); ++itr) pl->setEntry(itr->first, itr->second);
321  result = pl;
322  }
323  }
324 
325  return result;
326 }
327 
329 bool SIMPLEPreconditionerFactory::updateRequestedParameters(const Teuchos::ParameterList& pl) {
330  Teko_DEBUG_SCOPE("InvLSCStrategy::updateRequestedParameters", 10);
331  bool result = true;
332 
333  // update requested parameters in solvers
334  result &= invVelFactory_->updateRequestedParameters(pl);
335  result &= invPrsFactory_->updateRequestedParameters(pl);
336  if (customHFactory_ != Teuchos::null) result &= customHFactory_->updateRequestedParameters(pl);
337 
338  return result;
339 }
340 
341 } // end namespace NS
342 } // 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.
Preconditioner factory for building explcit inverse of diagonal operators. This includes block operat...
An implementation of a state object for block preconditioners.
virtual void setMassMatrix(Teko::LinearOp &mass)
Set the mass matrix for this factory.
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.
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.
virtual Teuchos::RCP< Teuchos::ParameterList > getRequestedParameters() const
For assisting in construction of the preconditioner.
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.