Teko  Version of the Day
 All Classes Files Functions Variables Pages
Teko_LSCPreconditionerFactory.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_LSCPreconditionerFactory.hpp"
48 
49 #include "Thyra_DefaultMultipliedLinearOp.hpp"
50 #include "Thyra_DefaultAddedLinearOp.hpp"
51 #include "Thyra_DefaultIdentityLinearOp.hpp"
52 #include "Thyra_DefaultZeroLinearOp.hpp"
53 #include "Thyra_get_Epetra_Operator.hpp"
54 
55 #include "Teko_LU2x2InverseOp.hpp"
56 #include "Teko_Utilities.hpp"
57 #include "Teko_BlockUpperTriInverseOp.hpp"
58 #include "Teko_StaticLSCStrategy.hpp"
59 #include "Teko_InvLSCStrategy.hpp"
60 #include "Teko_PresLaplaceLSCStrategy.hpp"
61 #include "Teko_LSCSIMPLECStrategy.hpp"
62 
63 #include "EpetraExt_RowMatrixOut.h"
64 
65 #include "Teuchos_Time.hpp"
66 
67 namespace Teko {
68 namespace NS {
69 
70 using Teuchos::rcp;
71 using Teuchos::rcp_dynamic_cast;
72 using Teuchos::RCP;
73 
74 using Thyra::multiply;
75 using Thyra::add;
76 using Thyra::identity;
77 
78 // Stabilized constructor
79 LSCPreconditionerFactory::LSCPreconditionerFactory(const LinearOp & invF,const LinearOp & invBQBtmC,
80  const LinearOp & invD,const LinearOp & invMass)
81  : invOpsStrategy_(rcp(new StaticLSCStrategy(invF,invBQBtmC,invD,invMass))), isSymmetric_(true)
82 { }
83 
84 // Stable constructor
85 LSCPreconditionerFactory::LSCPreconditionerFactory(const LinearOp & invF, const LinearOp & invBQBtmC,
86  const LinearOp & invMass)
87  : invOpsStrategy_(rcp(new StaticLSCStrategy(invF,invBQBtmC,invMass))), isSymmetric_(true)
88 { }
89 
90 // fully generic constructor
91 LSCPreconditionerFactory::LSCPreconditionerFactory(const RCP<LSCStrategy> & strategy)
92  : invOpsStrategy_(strategy), isSymmetric_(true)
93 { }
94 
95 LSCPreconditionerFactory::LSCPreconditionerFactory() : isSymmetric_(true)
96 { }
97 
98 // for PreconditionerFactoryBase
100 
101 // initialize a newly created preconditioner object
102 LinearOp LSCPreconditionerFactory::buildPreconditionerOperator(BlockedLinearOp & blockOp,BlockPreconditionerState & state) const
103 {
104  Teko_DEBUG_SCOPE("LSCPreconditionerFactory::buildPreconditionerOperator",10);
105  Teko_DEBUG_EXPR(Teuchos::Time timer(""));
106  Teko_DEBUG_EXPR(Teuchos::Time totalTimer(""));
107  Teko_DEBUG_EXPR(totalTimer.start());
108 
109  // extract sub-matrices from source operator
110  LinearOp F = blockOp->getBlock(0,0);
111  LinearOp B = blockOp->getBlock(1,0);
112  LinearOp Bt = blockOp->getBlock(0,1);
113 
114  if(not isSymmetric_)
115  Bt = scale(-1.0,adjoint(B));
116 
117  // build what is neccessary for the state object
118  Teko_DEBUG_EXPR(timer.start(true));
119  invOpsStrategy_->buildState(blockOp,state);
120  Teko_DEBUG_EXPR(timer.stop());
121  Teko_DEBUG_MSG("LSCPrecFact::buildPO BuildStateTime = " << timer.totalElapsedTime(),2);
122 
123  // extract operators from strategy
124  Teko_DEBUG_EXPR(timer.start(true));
125  LinearOp invF = invOpsStrategy_->getInvF(blockOp,state);
126  LinearOp invBQBtmC = invOpsStrategy_->getInvBQBt(blockOp,state);
127  LinearOp invBHBtmC = invOpsStrategy_->getInvBHBt(blockOp,state);
128  LinearOp outerStab = invOpsStrategy_->getOuterStabilization(blockOp,state);
129  LinearOp innerStab = invOpsStrategy_->getInnerStabilization(blockOp,state);
130 
131  // if necessary build an identity mass matrix
132  LinearOp invMass = invOpsStrategy_->getInvMass(blockOp,state);
133  LinearOp HScaling = invOpsStrategy_->getHScaling(blockOp,state);
134  if(invMass==Teuchos::null) invMass = identity<double>(F->range());
135  if(HScaling==Teuchos::null) HScaling = identity<double>(F->range());
136  Teko_DEBUG_EXPR(timer.stop());
137  Teko_DEBUG_MSG("LSCPrecFact::buildPO GetInvTime = " << timer.totalElapsedTime(),2);
138 
139  // need to build Schur complement, inv(P) = inv(B*Bt)*(B*F*Bt)*inv(B*Bt)
140 
141  // first construct middle operator: M = B * inv(Mass) * F * inv(Mass) * Bt
142  LinearOp M =
143  // (B * inv(Mass) ) * F * (inv(Mass) * Bt)
144  multiply( multiply(B,invMass), F , multiply(HScaling,Bt));
145  if(innerStab!=Teuchos::null) // deal with inner stabilization
146  M = add(M, innerStab);
147 
148  // now construct a linear operator schur complement
149  LinearOp invPschur;
150  if(outerStab!=Teuchos::null)
151  invPschur = add(multiply(invBQBtmC, M , invBHBtmC), outerStab);
152  else
153  invPschur = multiply(invBQBtmC, M , invBHBtmC);
154 
155  // build the preconditioner operator: Use LDU or upper triangular approximation
156  if(invOpsStrategy_->useFullLDU()) {
157  Teko_DEBUG_EXPR(totalTimer.stop());
158  Teko_DEBUG_MSG("LSCPrecFact::buildPO TotalTime = " << totalTimer.totalElapsedTime(),2);
159 
160  // solve using a full LDU decomposition
161  return createLU2x2InverseOp(blockOp,invF,invPschur,"LSC-LDU");
162  } else {
163  // build diagonal operations
164  std::vector<LinearOp> invDiag(2);
165  invDiag[0] = invF;
166  invDiag[1] = Thyra::scale(-1.0,invPschur);
167 
168  // get upper triangular matrix
169  BlockedLinearOp U = getUpperTriBlocks(blockOp);
170 
171  Teko_DEBUG_EXPR(totalTimer.stop());
172  Teko_DEBUG_MSG("LSCPrecFact::buildPO TotalTime = " << totalTimer.totalElapsedTime(),2);
173 
174  // solve using only one inversion of F
175  return createBlockUpperTriInverseOp(U,invDiag,"LSC-Upper");
176  }
177 }
178 
180 void LSCPreconditionerFactory::initializeFromParameterList(const Teuchos::ParameterList & pl)
181 {
182  Teko_DEBUG_SCOPE("LSCPreconditionerFactory::initializeFromParameterList",10);
183 
184  RCP<const InverseLibrary> invLib = getInverseLibrary();
185 
186  if(pl.isParameter("Is Symmetric"))
187  isSymmetric_ = pl.get<bool>("Is Symmetric");
188 
189  std::string name = "Basic Inverse";
190  if(pl.isParameter("Strategy Name"))
191  name = pl.get<std::string>("Strategy Name");
192  const Teuchos::ParameterEntry * pe = pl.getEntryPtr("Strategy Settings");
193 
194  // check for a mistake in input file
195  if(name!="Basic Inverse" && pe==0) {
196  RCP<Teuchos::FancyOStream> out = getOutputStream();
197  *out << "LSC Construction failed: ";
198  *out << "Strategy \"" << name << "\" requires a \"Strategy Settings\" sublist" << std::endl;
199  throw std::runtime_error("LSC Construction failed: Strategy Settings not set");
200  }
201 
202  // get the parameter list to construct the strategy
203  Teuchos::RCP<const Teuchos::ParameterList> stratPL = Teuchos::rcpFromRef(pl);
204  if(pe!=0)
205  stratPL = Teuchos::rcpFromRef(pl.sublist("Strategy Settings"));
206 
207  // build the strategy object
208  RCP<LSCStrategy> strategy = buildStrategy(name,*stratPL,invLib,getRequestHandler());
209 
210  // strategy could not be built
211  if(strategy==Teuchos::null) {
212  RCP<Teuchos::FancyOStream> out = getOutputStream();
213  *out << "LSC Construction failed: ";
214  *out << "Strategy \"" << name << "\" could not be constructed" << std::endl;
215  throw std::runtime_error("LSC Construction failed: Strategy could not be constructed");
216  }
217 
218  strategy->setSymmetric(isSymmetric_);
219  invOpsStrategy_ = strategy;
220 }
221 
223 Teuchos::RCP<Teuchos::ParameterList> LSCPreconditionerFactory::getRequestedParameters() const
224 {
225  return invOpsStrategy_->getRequestedParameters();
226 }
227 
229 bool LSCPreconditionerFactory::updateRequestedParameters(const Teuchos::ParameterList & pl)
230 {
231  return invOpsStrategy_->updateRequestedParameters(pl);
232 }
233 
235 // Static members and methods
237 
239 CloneFactory<LSCStrategy> LSCPreconditionerFactory::strategyBuilder_;
240 
253 RCP<LSCStrategy> LSCPreconditionerFactory::buildStrategy(const std::string & name,
254  const Teuchos::ParameterList & settings,
255  const RCP<const InverseLibrary> & invLib,
256  const RCP<RequestHandler> & rh)
257 {
258  Teko_DEBUG_SCOPE("LSCPreconditionerFactory::buildStrategy",10);
259 
260  // initialize the defaults if necessary
261  if(strategyBuilder_.cloneCount()==0) initializeStrategyBuilder();
262 
263  // request the preconditioner factory from the CloneFactory
264  Teko_DEBUG_MSG("Building LSC strategy \"" << name << "\"",1);
265  RCP<LSCStrategy> strategy = strategyBuilder_.build(name);
266 
267  if(strategy==Teuchos::null) return Teuchos::null;
268 
269  // now that inverse library has been set,
270  // pass in the parameter list
271  strategy->setRequestHandler(rh);
272  strategy->initializeFromParameterList(settings,*invLib);
273 
274  return strategy;
275 }
276 
290 void LSCPreconditionerFactory::addStrategy(const std::string & name,const RCP<Cloneable> & clone)
291 {
292  // initialize the defaults if necessary
293  if(strategyBuilder_.cloneCount()==0) initializeStrategyBuilder();
294 
295  // add clone to builder
296  strategyBuilder_.addClone(name,clone);
297 }
298 
300 void LSCPreconditionerFactory::initializeStrategyBuilder()
301 {
302  RCP<Cloneable> clone;
303 
304  // add various strategies to the factory
305  clone = rcp(new AutoClone<InvLSCStrategy>());
306  strategyBuilder_.addClone("Basic Inverse",clone);
307 
308  // add various strategies to the factory
309  clone = rcp(new AutoClone<PresLaplaceLSCStrategy>());
310  strategyBuilder_.addClone("Pressure Laplace",clone);
311 
312  // add various strategies to the factory
313  clone = rcp(new AutoClone<LSCSIMPLECStrategy>());
314  strategyBuilder_.addClone("SIMPLEC",clone);
315 }
316 
317 } // end namespace NS
318 } // end namespace Teko