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