47 #include "Teko_TimingsSIMPLEPreconditionerFactory.hpp" 50 #include "Teko_InverseFactory.hpp" 51 #include "Teko_BlockLowerTriInverseOp.hpp" 52 #include "Teko_BlockUpperTriInverseOp.hpp" 53 #include "Teko_DiagonalPreconditionerFactory.hpp" 55 #include "Teuchos_Time.hpp" 56 #include "Epetra_FECrsGraph.h" 57 #include "Epetra_FECrsMatrix.h" 58 #include "EpetraExt_PointToBlockDiagPermute.h" 59 #include "EpetraExt_MatrixMatrix.h" 60 #include "Thyra_EpetraOperatorWrapper.hpp" 61 #include "Thyra_EpetraLinearOp.hpp" 73 : invVelFactory_(inverse), invPrsFactory_(inverse), alpha_(alpha), fInverseType_(
Diagonal), useMass_(false)
74 , constrTotal_(
"SIMPLE Constr: Total"), subTotal_(
"SIMPLE Constr: Subs"), constrCount_(0)
79 const RCP<InverseFactory> & invPFact,
81 : invVelFactory_(invVFact), invPrsFactory_(invPFact), alpha_(alpha), fInverseType_(
Diagonal), useMass_(false)
82 , constrTotal_(
"SIMPLE Constr: Total"), subTotal_(
"SIMPLE Constr: Subs"), constrCount_(0)
86 : alpha_(1.0), fInverseType_(
Diagonal), useMass_(false)
87 , constrTotal_(
"SIMPLE Constr: Total"), subTotal_(
"SIMPLE Constr: Subs"), constrCount_(0)
92 if(constrTotal_.totalElapsedTime()>0.0) {
93 Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
94 out.setOutputToRootOnly(0);
96 out <<
"===========================================================================" << std::endl;
98 out <<
"SIMPLE Construction Count = " << constrCount_ << std::endl;
99 out <<
"SIMPLE Construction Total = " << constrTotal_.totalElapsedTime() << std::endl;
100 out <<
"SIMPLE Sub Components Total = " << subTotal_.totalElapsedTime() << std::endl;
102 out <<
"===========================================================================" << std::endl;
111 constrTotal_.start();
113 Teko_DEBUG_SCOPE(
"TimingsSIMPLEPreconditionerFactory::buildPreconditionerOperator",10);
114 Teko_DEBUG_EXPR(Teuchos::Time timer(
""));
119 TEUCHOS_ASSERT(rows==2);
120 TEUCHOS_ASSERT(cols==2);
122 bool buildExplicitSchurComplement =
true;
125 const LinearOp F =
getBlock(0,0,blockOp);
126 const LinearOp Bt =
getBlock(0,1,blockOp);
127 const LinearOp B =
getBlock(1,0,blockOp);
128 const LinearOp C =
getBlock(1,1,blockOp);
132 TEUCHOS_ASSERT(massMatrix_!=Teuchos::null);
137 std::string fApproxStr =
"<error>";
141 fApproxStr = customHFactory_->toString();
144 buildExplicitSchurComplement =
false;
146 else if(fInverseType_==
BlkDiag) {
153 buildExplicitSchurComplement =
true;
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);
virtual bool updateRequestedParameters(const Teuchos::ParameterList &pl)
For assisting in construction of the preconditioner.
virtual ~TimingsSIMPLEPreconditionerFactory()
Destructor that outputs construction timings.
Specifies that a block diagonal approximation is to be used.
Teuchos::RCP< RequestHandler > getRequestHandler() const
Get the request handler with pointers to the appropriate callbacks.
TimingsSIMPLEPreconditionerFactory()
Default constructor.
Specifies that just the diagonal is used.
int blockRowCount(const BlockedLinearOp &blo)
Get the row count in a block linear operator.
InverseLinearOp buildInverse(const InverseFactory &factory, const LinearOp &A)
Build an inverse operator using a factory and a linear operator.
virtual void setMassMatrix(Teko::LinearOp &mass)
Set the mass matrix for this factory.
Teuchos::RCP< const InverseLibrary > getInverseLibrary() const
Get the inverse library used by this preconditioner factory.
void scale(const double alpha, MultiVector &x)
Scale a multivector by a constant.
LinearOp buildPreconditionerOperator(BlockedLinearOp &blo, BlockPreconditionerState &state) const
Function that is called to build the preconditioner for the linear operator that is passed in...
LinearOp buildPreconditionerOperator(LinearOp &lo, PreconditionerState &state) const
void endBlockFill(BlockedLinearOp &blo)
Notify the blocked operator that the fill stage is completed.
virtual void initializeFromParameterList(const Teuchos::ParameterList &pl)
Initialize from a parameter list.
Preconditioner factory for building explcit inverse of diagonal operators. This includes block operat...
void rebuildInverse(const InverseFactory &factory, const LinearOp &A, InverseLinearOp &invA)
An implementation of a state object for block preconditioners.
VectorSpace rangeSpace(const LinearOp &lo)
Get the range space of a linear operator.
virtual Teuchos::RCP< Teuchos::ParameterList > getRequestedParameters() const
For assisting in construction of the preconditioner.
This linear operator prints diagnostics about operator application and creation times. It is useful for debugging problems and determining bottle necks.
MultiVector getBlock(int i, const BlockedMultiVector &bmv)
Get the ith block from a BlockedMultiVector object.
For user convenience, if Teko recieves this value, exceptions will be thrown.
virtual void initializeFromParameterList(const Teuchos::ParameterList &pl)
Initialize from a parameter list.
int blockColCount(const BlockedLinearOp &blo)
Get the column count in a block linear operator.
void setBlock(int i, int j, BlockedLinearOp &blo, const LinearOp &lo)
Set the i,j block in a BlockedLinearOp object.
virtual Teko::ModifiableLinearOp & getModifiableOp(const std::string &name)
Add a named operator to the state object.