47#include "Teko_TimingsSIMPLEPreconditionerFactory.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"
57#include "Teuchos_Time.hpp"
65TimingsSIMPLEPreconditionerFactory
66 ::TimingsSIMPLEPreconditionerFactory(
const RCP<InverseFactory> & inverse,
68 : invVelFactory_(inverse), invPrsFactory_(inverse), alpha_(alpha), fInverseType_(
Diagonal), useMass_(false)
69 , constrTotal_(
"SIMPLE Constr: Total"), subTotal_(
"SIMPLE Constr: Subs"), constrCount_(0)
72TimingsSIMPLEPreconditionerFactory
73 ::TimingsSIMPLEPreconditionerFactory(
const RCP<InverseFactory> & invVFact,
74 const RCP<InverseFactory> & invPFact,
76 : invVelFactory_(invVFact), invPrsFactory_(invPFact), alpha_(alpha), fInverseType_(
Diagonal), useMass_(false)
77 , constrTotal_(
"SIMPLE Constr: Total"), subTotal_(
"SIMPLE Constr: Subs"), constrCount_(0)
80TimingsSIMPLEPreconditionerFactory::TimingsSIMPLEPreconditionerFactory()
81 : alpha_(1.0), fInverseType_(
Diagonal), useMass_(false)
82 , constrTotal_(
"SIMPLE Constr: Total"), subTotal_(
"SIMPLE Constr: Subs"), constrCount_(0)
87 if(constrTotal_.totalElapsedTime()>0.0) {
88 Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
89 out.setOutputToRootOnly(0);
91 out <<
"===========================================================================" << std::endl;
93 out <<
"SIMPLE Construction Count = " << constrCount_ << std::endl;
94 out <<
"SIMPLE Construction Total = " << constrTotal_.totalElapsedTime() << std::endl;
95 out <<
"SIMPLE Sub Components Total = " << subTotal_.totalElapsedTime() << std::endl;
97 out <<
"===========================================================================" << std::endl;
102LinearOp TimingsSIMPLEPreconditionerFactory
103 ::buildPreconditionerOperator(BlockedLinearOp & blockOp,
106 constrTotal_.start();
108 Teko_DEBUG_SCOPE(
"TimingsSIMPLEPreconditionerFactory::buildPreconditionerOperator",10);
109 Teko_DEBUG_EXPR(Teuchos::Time timer(
""));
114 TEUCHOS_ASSERT(rows==2);
115 TEUCHOS_ASSERT(cols==2);
117 bool buildExplicitSchurComplement =
true;
120 const LinearOp F =
getBlock(0,0,blockOp);
121 const LinearOp Bt =
getBlock(0,1,blockOp);
122 const LinearOp B =
getBlock(1,0,blockOp);
123 const LinearOp C =
getBlock(1,1,blockOp);
127 TEUCHOS_ASSERT(massMatrix_!=Teuchos::null);
132 std::string fApproxStr =
"<error>";
136 fApproxStr = customHFactory_->toString();
139 buildExplicitSchurComplement =
false;
141 else if(fInverseType_==
BlkDiag) {
142#ifdef TEKO_HAVE_EPETRA
149 buildExplicitSchurComplement =
true;
152 throw std::logic_error(
"TimingsSIMPLEPreconditionerFactory fInverseType_ == "
153 "BlkDiag but EPETRA is turned off!");
159 H = getInvDiagonalOp(matF,fInverseType_);
161 fApproxStr = getDiagonalName(fInverseType_);
166 RCP<const Teuchos::ParameterList> pl = state.getParameterList();
168 if(pl->isParameter(
"stepsize")) {
170 double stepsize = pl->get<
double>(
"stepsize");
174 H =
scale(stepsize,H);
181 if(buildExplicitSchurComplement) {
188 mHBt = explicitMultiply(H,Bt,mHBt);
194 BHBt = explicitMultiply(B,HBt,BHBt);
199 mhatS = explicitAdd(C,
scale(-1.0,BHBt),mhatS);
205 HBt = multiply(H,Bt);
207 hatS = add(C,
scale(-1.0,multiply(B,HBt)));
211 if(timed_HBt_==Teuchos::null) {
215 timed_HBt_->setLinearOp(HBt);
219 if(timed_B_==Teuchos::null) {
223 timed_B_->setLinearOp(B);
229 if(invF==Teuchos::null) {
237 timed_invF_->setLinearOp(invF);
244 if(invS==Teuchos::null) {
252 timed_invS_->setLinearOp(invS);
256 std::vector<LinearOp> invDiag(2);
259 BlockedLinearOp L = zeroBlockedOp(blockOp);
263 invDiag[0] = timed_invF_;
264 invDiag[1] = timed_invS_;
265 LinearOp invL = createBlockLowerTriInverseOp(L,invDiag);
268 BlockedLinearOp U = zeroBlockedOp(blockOp);
269 setBlock(0,1,U,scale<double>(1.0/alpha_,timed_HBt_.getConst()));
274 LinearOp invU = createBlockUpperTriInverseOp(U,invDiag);
277 Teko::LinearOp iU_t_iL = multiply(invU,invL,
"SIMPLE_"+fApproxStr);
280 if(timed_iU_t_iL_==Teuchos::null)
281 timed_iU_t_iL_ = Teuchos::rcp(
new DiagnosticLinearOp(getOutputStream(),iU_t_iL,
"iU_t_iL"));
283 timed_iU_t_iL_->setLinearOp(iU_t_iL);
289 return timed_iU_t_iL_;
299 customHFactory_ = Teuchos::null;
303 std::string invStr=
"", invVStr=
"", invPStr=
"";
307 if(pl.isParameter(
"Inverse Type"))
308 invStr = pl.get<std::string>(
"Inverse Type");
309 if(pl.isParameter(
"Inverse Velocity Type"))
310 invVStr = pl.get<std::string>(
"Inverse Velocity Type");
311 if(pl.isParameter(
"Inverse Pressure Type"))
312 invPStr = pl.get<std::string>(
"Inverse Pressure Type");
313 if(pl.isParameter(
"Alpha"))
314 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");
319 fInverseType_ = getDiagonalType(fInverseStr);
321 customHFactory_ = invLib->getInverseFactory(fInverseStr);
325 BlkDiagList_=pl.sublist(
"H options");
327 if(pl.isParameter(
"Use Mass Scaling"))
328 useMass_ = pl.get<
bool>(
"Use Mass Scaling");
330 Teko_DEBUG_MSG_BEGIN(5)
331 DEBUG_STREAM <<
"SIMPLE Parameters: " << std::endl;
332 DEBUG_STREAM <<
" inv type = \"" << invStr <<
"\"" << std::endl;
333 DEBUG_STREAM <<
" inv v type = \"" << invVStr <<
"\"" << std::endl;
334 DEBUG_STREAM <<
" inv p type = \"" << invPStr <<
"\"" << std::endl;
335 DEBUG_STREAM <<
" alpha = " << alpha_ << std::endl;
336 DEBUG_STREAM <<
" use mass = " << useMass_ << std::endl;
337 DEBUG_STREAM <<
" vel scaling = " << getDiagonalName(fInverseType_) << std::endl;
338 DEBUG_STREAM <<
"SIMPLE Parameter list: " << std::endl;
339 pl.print(DEBUG_STREAM);
343 if(invStr==
"") invStr =
"Amesos";
344 if(invVStr==
"") invVStr = invStr;
345 if(invPStr==
"") invPStr = invStr;
348 RCP<InverseFactory> invVFact, invPFact;
351 invVFact = invLib->getInverseFactory(invVStr);
354 invPFact = invLib->getInverseFactory(invPStr);
357 invVelFactory_ = invVFact;
358 invPrsFactory_ = invPFact;
362 rh->preRequest<Teko::LinearOp>(Teko::RequestMesg(
"Velocity Mass Matrix"));
364 = rh->request<Teko::LinearOp>(Teko::RequestMesg(
"Velocity Mass Matrix"));
372 Teuchos::RCP<Teuchos::ParameterList> result;
373 Teuchos::RCP<Teuchos::ParameterList> pl = rcp(
new Teuchos::ParameterList());
376 RCP<Teuchos::ParameterList> vList = invVelFactory_->getRequestedParameters();
377 if(vList!=Teuchos::null) {
378 Teuchos::ParameterList::ConstIterator itr;
379 for(itr=vList->begin();itr!=vList->end();++itr)
380 pl->setEntry(itr->first,itr->second);
385 RCP<Teuchos::ParameterList> pList = invPrsFactory_->getRequestedParameters();
386 if(pList!=Teuchos::null) {
387 Teuchos::ParameterList::ConstIterator itr;
388 for(itr=pList->begin();itr!=pList->end();++itr)
389 pl->setEntry(itr->first,itr->second);
394 if(customHFactory_!=Teuchos::null) {
395 RCP<Teuchos::ParameterList> hList = customHFactory_->getRequestedParameters();
396 if(hList!=Teuchos::null) {
397 Teuchos::ParameterList::ConstIterator itr;
398 for(itr=hList->begin();itr!=hList->end();++itr)
399 pl->setEntry(itr->first,itr->second);
410 Teko_DEBUG_SCOPE(
"InvLSCStrategy::updateRequestedParameters",10);
414 result &= invVelFactory_->updateRequestedParameters(pl);
415 result &= invPrsFactory_->updateRequestedParameters(pl);
416 if(customHFactory_!=Teuchos::null)
417 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.
This linear operator prints diagnostics about operator application and creation times....
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 ~TimingsSIMPLEPreconditionerFactory()
Destructor that outputs construction timings.
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.