47#include "Teko_SIMPLEPreconditionerFactory.hpp"
50#include "Teko_InverseFactory.hpp"
51#include "Teko_BlockLowerTriInverseOp.hpp"
52#include "Teko_BlockUpperTriInverseOp.hpp"
54#ifdef TEKO_HAVE_EPETRA
55#include "Teko_DiagonalPreconditionerFactory.hpp"
58#include "Teuchos_Time.hpp"
66SIMPLEPreconditionerFactory
67 ::SIMPLEPreconditionerFactory(
const RCP<InverseFactory> & inverse,
69 : invVelFactory_(inverse), invPrsFactory_(inverse), alpha_(alpha), fInverseType_(
Diagonal), useMass_(false)
72SIMPLEPreconditionerFactory
73 ::SIMPLEPreconditionerFactory(
const RCP<InverseFactory> & invVFact,
74 const RCP<InverseFactory> & invPFact,
76 : invVelFactory_(invVFact), invPrsFactory_(invPFact), alpha_(alpha), fInverseType_(
Diagonal), useMass_(false)
79SIMPLEPreconditionerFactory::SIMPLEPreconditionerFactory()
80 : alpha_(1.0), fInverseType_(
Diagonal), useMass_(false)
84LinearOp SIMPLEPreconditionerFactory
85 ::buildPreconditionerOperator(BlockedLinearOp & blockOp,
88 Teko_DEBUG_SCOPE(
"SIMPLEPreconditionerFactory::buildPreconditionerOperator",10);
89 Teko_DEBUG_EXPR(Teuchos::Time timer(
""));
94 TEUCHOS_ASSERT(rows==2);
95 TEUCHOS_ASSERT(cols==2);
97 bool buildExplicitSchurComplement =
true;
100 const LinearOp F =
getBlock(0,0,blockOp);
101 const LinearOp Bt =
getBlock(0,1,blockOp);
102 const LinearOp B =
getBlock(1,0,blockOp);
103 const LinearOp C =
getBlock(1,1,blockOp);
107 TEUCHOS_ASSERT(massMatrix_!=Teuchos::null);
112 std::string fApproxStr =
"<error>";
116 fApproxStr = customHFactory_->toString();
119 buildExplicitSchurComplement =
false;
121 else if(fInverseType_==
BlkDiag) {
122#ifdef TEKO_HAVE_EPETRA
135 buildExplicitSchurComplement =
true;
138 throw std::logic_error(
"SIMPLEPreconditionerFactory fInverseType_ == "
139 "BlkDiag but EPETRA is turned off!");
145 H = getInvDiagonalOp(matF,fInverseType_);
146 fApproxStr = getDiagonalName(fInverseType_);
151 RCP<const Teuchos::ParameterList> pl = state.getParameterList();
153 if(pl->isParameter(
"stepsize")) {
155 double stepsize = pl->get<
double>(
"stepsize");
159 H =
scale(stepsize,H);
166 if(buildExplicitSchurComplement) {
172 mHBt = explicitMultiply(H,Bt,mHBt);
176 BHBt = explicitMultiply(B,HBt,BHBt);
179 mhatS = explicitAdd(C,
scale(-1.0,BHBt),mhatS);
184 HBt = multiply(H,Bt);
186 hatS = add(C,
scale(-1.0,multiply(B,HBt)));
191 if(invF==Teuchos::null)
198 if(invS==Teuchos::null)
203 std::vector<LinearOp> invDiag(2);
206 BlockedLinearOp L = zeroBlockedOp(blockOp);
212 LinearOp invL = createBlockLowerTriInverseOp(L,invDiag);
215 BlockedLinearOp U = zeroBlockedOp(blockOp);
221 LinearOp invU = createBlockUpperTriInverseOp(U,invDiag);
224 return multiply(invU,invL,
"SIMPLE_"+fApproxStr);
234 customHFactory_ = Teuchos::null;
238 std::string invStr=
"", invVStr=
"", invPStr=
"";
242 if(pl.isParameter(
"Inverse Type"))
243 invStr = pl.get<std::string>(
"Inverse Type");
244 if(pl.isParameter(
"Inverse Velocity Type"))
245 invVStr = pl.get<std::string>(
"Inverse Velocity Type");
246 if(pl.isParameter(
"Inverse Pressure Type"))
247 invPStr = pl.get<std::string>(
"Inverse Pressure Type");
248 if(pl.isParameter(
"Alpha"))
249 alpha_ = pl.get<
double>(
"Alpha");
250 if(pl.isParameter(
"Explicit Velocity Inverse Type")) {
251 std::string fInverseStr = pl.get<std::string>(
"Explicit Velocity Inverse Type");
254 fInverseType_ = getDiagonalType(fInverseStr);
256 customHFactory_ = invLib->getInverseFactory(fInverseStr);
260 BlkDiagList_=pl.sublist(
"H options");
262 if(pl.isParameter(
"Use Mass Scaling"))
263 useMass_ = pl.get<
bool>(
"Use Mass Scaling");
265 Teko_DEBUG_MSG_BEGIN(5)
266 DEBUG_STREAM <<
"SIMPLE Parameters: " << std::endl;
267 DEBUG_STREAM <<
" inv type = \"" << invStr <<
"\"" << std::endl;
268 DEBUG_STREAM <<
" inv v type = \"" << invVStr <<
"\"" << std::endl;
269 DEBUG_STREAM <<
" inv p type = \"" << invPStr <<
"\"" << std::endl;
270 DEBUG_STREAM <<
" alpha = " << alpha_ << std::endl;
271 DEBUG_STREAM <<
" use mass = " << useMass_ << std::endl;
272 DEBUG_STREAM <<
" vel scaling = " << getDiagonalName(fInverseType_) << std::endl;
273 DEBUG_STREAM <<
"SIMPLE Parameter list: " << std::endl;
274 pl.print(DEBUG_STREAM);
278 if(invStr==
"") invStr =
"Amesos";
279 if(invVStr==
"") invVStr = invStr;
280 if(invPStr==
"") invPStr = invStr;
283 RCP<InverseFactory> invVFact, invPFact;
286 invVFact = invLib->getInverseFactory(invVStr);
289 invPFact = invLib->getInverseFactory(invPStr);
292 invVelFactory_ = invVFact;
293 invPrsFactory_ = invPFact;
297 rh->preRequest<Teko::LinearOp>(Teko::RequestMesg(
"Velocity Mass Matrix"));
299 = rh->request<Teko::LinearOp>(Teko::RequestMesg(
"Velocity Mass Matrix"));
307 Teuchos::RCP<Teuchos::ParameterList> result;
308 Teuchos::RCP<Teuchos::ParameterList> pl = rcp(
new Teuchos::ParameterList());
311 RCP<Teuchos::ParameterList> vList = invVelFactory_->getRequestedParameters();
312 if(vList!=Teuchos::null) {
313 Teuchos::ParameterList::ConstIterator itr;
314 for(itr=vList->begin();itr!=vList->end();++itr)
315 pl->setEntry(itr->first,itr->second);
320 RCP<Teuchos::ParameterList> pList = invPrsFactory_->getRequestedParameters();
321 if(pList!=Teuchos::null) {
322 Teuchos::ParameterList::ConstIterator itr;
323 for(itr=pList->begin();itr!=pList->end();++itr)
324 pl->setEntry(itr->first,itr->second);
329 if(customHFactory_!=Teuchos::null) {
330 RCP<Teuchos::ParameterList> hList = customHFactory_->getRequestedParameters();
331 if(hList!=Teuchos::null) {
332 Teuchos::ParameterList::ConstIterator itr;
333 for(itr=hList->begin();itr!=hList->end();++itr)
334 pl->setEntry(itr->first,itr->second);
345 Teko_DEBUG_SCOPE(
"InvLSCStrategy::updateRequestedParameters",10);
349 result &= invVelFactory_->updateRequestedParameters(pl);
350 result &= invPrsFactory_->updateRequestedParameters(pl);
351 if(customHFactory_!=Teuchos::null)
352 result &= customHFactory_->updateRequestedParameters(pl);
InverseLinearOp buildInverse(const InverseFactory &factory, const LinearOp &A)
Build an inverse operator using a factory and a linear operator.
void rebuildInverse(const InverseFactory &factory, const LinearOp &A, InverseLinearOp &invA)
void scale(const double alpha, MultiVector &x)
Scale a multivector by a constant.
@ BlkDiag
Specifies that a block diagonal approximation is to be used.
@ NotDiag
For user convenience, if Teko recieves this value, exceptions will be thrown.
@ Diagonal
Specifies that just the diagonal is used.
int blockRowCount(const BlockedLinearOp &blo)
Get the row count in a block linear operator.
MultiVector getBlock(int i, const BlockedMultiVector &bmv)
Get the ith block from a BlockedMultiVector object.
void setBlock(int i, int j, BlockedLinearOp &blo, const LinearOp &lo)
Set the i,j block in a BlockedLinearOp object.
VectorSpace rangeSpace(const LinearOp &lo)
Replace nonzeros with a scalar value, used to zero out an operator.
int blockColCount(const BlockedLinearOp &blo)
Get the column count in a block linear operator.
void endBlockFill(BlockedLinearOp &blo)
Notify the blocked operator that the fill stage is completed.
An implementation of a state object for block preconditioners.
Preconditioner factory for building explcit inverse of diagonal operators. This includes block operat...
virtual void initializeFromParameterList(const Teuchos::ParameterList &pl)
Initialize from a parameter list.
LinearOp buildPreconditionerOperator(LinearOp &lo, PreconditionerState &state) const
virtual bool updateRequestedParameters(const Teuchos::ParameterList &pl)
For assisting in construction of the preconditioner.
virtual void initializeFromParameterList(const Teuchos::ParameterList &pl)
Initialize from a parameter list.
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.
Teuchos::RCP< const InverseLibrary > getInverseLibrary() const
Get the inverse library used by this preconditioner factory.
Teuchos::RCP< RequestHandler > getRequestHandler() const
Get the request handler with pointers to the appropriate callbacks.
virtual Teko::ModifiableLinearOp & getModifiableOp(const std::string &name)
Add a named operator to the state object.