47#include "NS/Teko_InvLSCStrategy.hpp"
49#include "Thyra_DefaultDiagonalLinearOp.hpp"
50#include "Thyra_VectorStdOps.hpp"
52#include "Teko_ConfigDefs.hpp"
54#ifdef TEKO_HAVE_EPETRA
55#include "Thyra_EpetraThyraWrappers.hpp"
56#include "Thyra_get_Epetra_Operator.hpp"
57#include "Thyra_EpetraLinearOp.hpp"
58#include "Epetra_Vector.h"
59#include "Epetra_Map.h"
60#include "EpetraExt_RowMatrixOut.h"
61#include "EpetraExt_MultiVectorOut.h"
62#include "EpetraExt_VectorOut.h"
63#include "Teko_EpetraHelpers.hpp"
64#include "Teko_EpetraOperatorWrapper.hpp"
67#include "Teuchos_Time.hpp"
68#include "Teuchos_TimeMonitor.hpp"
72#include "NS/Teko_LSCPreconditionerFactory.hpp"
74#include "Teko_TpetraHelpers.hpp"
76#include "Thyra_TpetraLinearOp.hpp"
79using Teuchos::rcp_dynamic_cast;
80using Teuchos::rcp_const_cast;
91InvLSCStrategy::InvLSCStrategy()
92 : massMatrix_(Teuchos::null), invFactoryF_(Teuchos::null), invFactoryS_(Teuchos::null), eigSolveParam_(5)
93 , rowZeroingNeeded_(false), useFullLDU_(false), useMass_(false), useLumping_(false), useWScaling_(false), scaleType_(
Diagonal)
94 , isSymmetric_(true), assumeStable_(false)
97InvLSCStrategy::InvLSCStrategy(
const Teuchos::RCP<InverseFactory> & factory,
bool rzn)
98 : massMatrix_(Teuchos::null), invFactoryF_(factory), invFactoryS_(factory), eigSolveParam_(5), rowZeroingNeeded_(rzn)
99 , useFullLDU_(false), useMass_(false), useLumping_(false), useWScaling_(false), scaleType_(
Diagonal)
100 , isSymmetric_(true), assumeStable_(false)
103InvLSCStrategy::InvLSCStrategy(
const Teuchos::RCP<InverseFactory> & invFactF,
104 const Teuchos::RCP<InverseFactory> & invFactS,
106 : massMatrix_(Teuchos::null), invFactoryF_(invFactF), invFactoryS_(invFactS), eigSolveParam_(5), rowZeroingNeeded_(rzn)
107 , useFullLDU_(false), useMass_(false), useLumping_(false), useWScaling_(false), scaleType_(
Diagonal)
108 , isSymmetric_(true), assumeStable_(false)
111InvLSCStrategy::InvLSCStrategy(
const Teuchos::RCP<InverseFactory> & factory,LinearOp & mass,
bool rzn)
112 : massMatrix_(mass), invFactoryF_(factory), invFactoryS_(factory), eigSolveParam_(5), rowZeroingNeeded_(rzn)
113 , useFullLDU_(false), useMass_(false), useLumping_(false), useWScaling_(false), scaleType_(
Diagonal)
114 , isSymmetric_(true), assumeStable_(false)
117InvLSCStrategy::InvLSCStrategy(
const Teuchos::RCP<InverseFactory> & invFactF,
118 const Teuchos::RCP<InverseFactory> & invFactS,
119 LinearOp & mass,
bool rzn)
120 : massMatrix_(mass), invFactoryF_(invFactF), invFactoryS_(invFactS), eigSolveParam_(5), rowZeroingNeeded_(rzn)
121 , useFullLDU_(false), useMass_(false), useLumping_(false), useWScaling_(false), scaleType_(
Diagonal)
122 , isSymmetric_(true), assumeStable_(false)
129 Teko_DEBUG_SCOPE(
"InvLSCStrategy::buildState",10);
132 TEUCHOS_ASSERT(lscState!=0);
136 Teko_DEBUG_EXPR(Teuchos::Time timer(
""));
140 Teko_DEBUG_SCOPE(
"LSC::buildState constructing operators",1);
141 Teko_DEBUG_EXPR(timer.start(
true));
143 initializeState(A,lscState);
145 Teko_DEBUG_EXPR(timer.stop());
146 Teko_DEBUG_MSG(
"LSC::buildState BuildOpsTime = " << timer.totalElapsedTime(),1);
151 Teko_DEBUG_SCOPE(
"LSC::buildState calculating inverses",1);
152 Teko_DEBUG_EXPR(timer.start(
true));
154 computeInverses(A,lscState);
156 Teko_DEBUG_EXPR(timer.stop());
157 Teko_DEBUG_MSG(
"LSC::buildState BuildInvTime = " << timer.totalElapsedTime(),1);
181 TEUCHOS_ASSERT(lscState!=0);
184 return lscState->
aiD_;
190 TEUCHOS_ASSERT(lscState!=0);
198 if(hScaling_!=Teuchos::null)
return hScaling_;
199 return getInvMass(A,state);
203void InvLSCStrategy::initializeState(
const BlockedLinearOp & A,
LSCPrecondState * state)
const
205 Teko_DEBUG_SCOPE(
"InvLSCStrategy::initializeState",10);
208 const LinearOp Bt =
getBlock(0,1,A);
213 LinearOp G = isSymmetric_ ? Bt :
adjoint(D);
215 bool isStabilized = assumeStable_ ? false : (not isZeroOp(C));
223 if(massMatrix_==Teuchos::null) {
224 Teko_DEBUG_MSG(
"LSC::initializeState Build Scaling <F> type \""
225 << getDiagonalName(scaleType_) <<
"\"" ,1);
226 state->
invMass_ = getInvDiagonalOp(F,scaleType_);
228 else if(state->
invMass_==Teuchos::null) {
229 Teko_DEBUG_MSG(
"LSC::initializeState Build Scaling <mass> type \""
230 << getDiagonalName(scaleType_) <<
"\"" ,1);
231 state->
invMass_ = getInvDiagonalOp(massMatrix_,scaleType_);
237 Teko_DEBUG_MSG(
"Computed BQBt",10);
240 if(wScaling_!=Teuchos::null && hScaling_==Teuchos::null) {
242 RCP<const Thyra::VectorBase<double> > w = wScaling_->col(0);
243 RCP<const Thyra::VectorBase<double> > iQu
244 = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<double> >(state->
invMass_)->getDiag();
245 RCP<Thyra::VectorBase<double> > h = Thyra::createMember(iQu->space());
247 Thyra::put_scalar(0.0,h.ptr());
248 Thyra::ele_wise_prod(1.0,*w,*iQu,h.ptr());
249 hScaling_ = Teuchos::rcp(
new Thyra::DefaultDiagonalLinearOp<double>(h));
252 LinearOp H = hScaling_;
253 if(H==Teuchos::null && not isSymmetric_)
260 RCP<Teuchos::Time> time = Teuchos::TimeMonitor::getNewTimer(
"InvLSCStrategy::initializeState Build BHBt");
261 Teuchos::TimeMonitor timer(*time);
264 state->
BHBt_ = explicitMultiply(D,H,G,state->
BHBt_);
268 if(not isStabilized) {
273 state->
aiD_ = Teuchos::null;
282 if(!Teko::TpetraHelpers::isTpetraLinearOp(F)){
283#ifdef TEKO_HAVE_EPETRA
284 const RCP<const Epetra_Operator> epF = Thyra::get_Epetra_Operator(*F);
285 if(epF!=Teuchos::null && rowZeroingNeeded_) {
287 const RCP<const Epetra_CrsMatrix> crsF = rcp_dynamic_cast<const Epetra_CrsMatrix>(epF);
290 if(crsF!=Teuchos::null) {
291 std::vector<int> zeroIndices;
294 Teko::Epetra::identityRowIndices(crsF->RowMap(), *crsF,zeroIndices);
301 throw std::logic_error(
"InvLSCStrategy::initializeState is trying to use "
302 "Epetra code, but TEKO is not built with Epetra!");
307 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > crsF = Teko::TpetraHelpers::getTpetraCrsMatrix(F, &scalar, &transp);
309 std::vector<GO> zeroIndices;
312 Teko::TpetraHelpers::identityRowIndices(*crsF->getRowMap(), *crsF,zeroIndices);
315 modF = Thyra::tpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(crsF->getDomainMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(crsF->getRangeMap()),rcp(
new Teko::TpetraHelpers::ZeroedOperator(zeroIndices,crsF)));
319 Teko_DEBUG_MSG(
"Calculating gamma",10);
320 LinearOp iQuF = multiply(state->
invMass_,modF);
323 Teko::LinearOp stabMatrix;
324 state->
gamma_ = std::fabs(Teko::computeSpectralRad(iQuF,5e-2,
false,eigSolveParam_))/3.0;
325 Teko_DEBUG_MSG(
"Calculated gamma",10);
326 if(userPresStabMat_!=Teuchos::null) {
327 Teko::LinearOp invDGl = Teko::getInvDiagonalOp(userPresStabMat_);
328 Teko::LinearOp gammaOp = multiply(invDGl,C);
329 state->
gamma_ *= std::fabs(Teko::computeSpectralRad(gammaOp,5e-2,
false,eigSolveParam_));
330 stabMatrix = userPresStabMat_;
336 LinearOp invDiagF = getInvDiagonalOp(F);
337 Teko::ModifiableLinearOp modB_idF_Bt = state->
getInverse(
"BidFBt");
338 modB_idF_Bt = explicitMultiply(B,invDiagF,Bt,modB_idF_Bt);
340 const LinearOp B_idF_Bt = modB_idF_Bt;
342 MultiVector vec_D = getDiagonal(B_idF_Bt);
343 update(-1.0,getDiagonal(C),1.0,vec_D);
344 const LinearOp invD = buildInvDiagonal(vec_D,
"inv(D)");
346 Teko_DEBUG_MSG(
"Calculating alpha",10);
347 const LinearOp BidFBtidD = multiply<double>(B_idF_Bt,invD);
348 double num = std::fabs(Teko::computeSpectralRad(BidFBtidD,5e-2,
false,eigSolveParam_));
349 Teko_DEBUG_MSG(
"Calculated alpha",10);
351 state->
aiD_ = Thyra::scale(state->
alpha_,invD);
354 Teko::ModifiableLinearOp BQBtmC = state->
getInverse(
"BQBtmC");
359 Teko::ModifiableLinearOp BHBtmC = state->
getInverse(
"BHBtmC");
367 Teko_DEBUG_MSG_BEGIN(5)
368 DEBUG_STREAM <<
"LSC Gamma Parameter = " << state->
gamma_ << std::endl;
369 DEBUG_STREAM <<
"LSC Alpha Parameter = " << state->
alpha_ << std::endl;
380void InvLSCStrategy::computeInverses(
const BlockedLinearOp & A,
LSCPrecondState * state)
const
382 Teko_DEBUG_SCOPE(
"InvLSCStrategy::computeInverses",10);
383 Teko_DEBUG_EXPR(Teuchos::Time invTimer(
""));
390 Teko_DEBUG_MSG(
"LSC::computeInverses Building inv(F)",1);
391 Teko_DEBUG_EXPR(invTimer.start(
true));
392 InverseLinearOp invF = state->
getInverse(
"invF");
393 if(invF==Teuchos::null) {
399 Teko_DEBUG_EXPR(invTimer.stop());
400 Teko_DEBUG_MSG(
"LSC::computeInverses GetInvF = " << invTimer.totalElapsedTime(),1);
405 Teko_DEBUG_MSG(
"LSC::computeInverses Building inv(BQBtmC)",1);
406 Teko_DEBUG_EXPR(invTimer.start(
true));
407 const LinearOp BQBt = state->
getInverse(
"BQBtmC");
408 InverseLinearOp invBQBt = state->
getInverse(
"invBQBtmC");
409 if(invBQBt==Teuchos::null) {
415 Teko_DEBUG_EXPR(invTimer.stop());
416 Teko_DEBUG_MSG(
"LSC::computeInverses GetInvBQBt = " << invTimer.totalElapsedTime(),1);
421 ModifiableLinearOp invBHBt = state->
getInverse(
"invBHBtmC");
422 if(hScaling_!=Teuchos::null || not isSymmetric_) {
424 Teko_DEBUG_MSG(
"LSC::computeInverses Building inv(BHBtmC)",1);
425 Teko_DEBUG_EXPR(invTimer.start(
true));
426 const LinearOp BHBt = state->
getInverse(
"BHBtmC");
427 if(invBHBt==Teuchos::null) {
433 Teko_DEBUG_EXPR(invTimer.stop());
434 Teko_DEBUG_MSG(
"LSC::computeInverses GetInvBHBt = " << invTimer.totalElapsedTime(),1);
436 else if(invBHBt==Teuchos::null) {
443void InvLSCStrategy::initializeFromParameterList(
const Teuchos::ParameterList & pl,
const InverseLibrary & invLib)
446 std::string invStr=
"", invVStr=
"", invPStr=
"";
447 bool rowZeroing =
true;
452 if(pl.isParameter(
"Inverse Type"))
453 invStr = pl.get<std::string>(
"Inverse Type");
454 if(pl.isParameter(
"Inverse Velocity Type"))
455 invVStr = pl.get<std::string>(
"Inverse Velocity Type");
456 if(pl.isParameter(
"Inverse Pressure Type"))
457 invPStr = pl.get<std::string>(
"Inverse Pressure Type");
458 if(pl.isParameter(
"Ignore Boundary Rows"))
459 rowZeroing = pl.get<
bool>(
"Ignore Boundary Rows");
460 if(pl.isParameter(
"Use LDU"))
461 useLDU = pl.get<
bool>(
"Use LDU");
462 if(pl.isParameter(
"Use Mass Scaling"))
463 useMass_ = pl.get<
bool>(
"Use Mass Scaling");
466 if(pl.isParameter(
"Use W-Scaling"))
467 useWScaling_ = pl.get<
bool>(
"Use W-Scaling");
468 if(pl.isParameter(
"Eigen Solver Iterations"))
469 eigSolveParam_ = pl.get<
int>(
"Eigen Solver Iterations");
470 if(pl.isParameter(
"Scaling Type")) {
471 scaleType_ = getDiagonalType(pl.get<std::string>(
"Scaling Type"));
472 TEUCHOS_TEST_FOR_EXCEPT(scaleType_==
NotDiag);
474 if(pl.isParameter(
"Assume Stable Discretization"))
475 assumeStable_ = pl.get<
bool>(
"Assume Stable Discretization");
477 Teko_DEBUG_MSG_BEGIN(5)
478 DEBUG_STREAM <<
"LSC Inverse Strategy Parameters: " << std::endl;
479 DEBUG_STREAM <<
" inv type = \"" << invStr <<
"\"" << std::endl;
480 DEBUG_STREAM <<
" inv v type = \"" << invVStr <<
"\"" << std::endl;
481 DEBUG_STREAM <<
" inv p type = \"" << invPStr <<
"\"" << std::endl;
482 DEBUG_STREAM <<
" bndry rows = " << rowZeroing << std::endl;
483 DEBUG_STREAM <<
" use ldu = " << useLDU << std::endl;
484 DEBUG_STREAM <<
" use mass = " << useMass_ << std::endl;
485 DEBUG_STREAM <<
" use w-scaling = " << useWScaling_ << std::endl;
486 DEBUG_STREAM <<
" assume stable = " << assumeStable_ << std::endl;
487 DEBUG_STREAM <<
" scale type = " << getDiagonalName(scaleType_) << std::endl;
488 DEBUG_STREAM <<
"LSC Inverse Strategy Parameter list: " << std::endl;
489 pl.print(DEBUG_STREAM);
493 if(invStr==
"") invStr =
"Amesos";
494 if(invVStr==
"") invVStr = invStr;
495 if(invPStr==
"") invPStr = invStr;
498 invFactoryF_ = invLib.getInverseFactory(invVStr);
499 invFactoryS_ = invFactoryF_;
501 invFactoryS_ = invLib.getInverseFactory(invPStr);
504 setUseFullLDU(useLDU);
505 setRowZeroing(rowZeroing);
508 Teuchos::RCP<Teko::RequestHandler> rh = getRequestHandler();
509 rh->preRequest<Teko::LinearOp>(Teko::RequestMesg(
"Velocity Mass Matrix"));
511 = rh->request<Teko::LinearOp>(Teko::RequestMesg(
"Velocity Mass Matrix"));
518Teuchos::RCP<Teuchos::ParameterList> InvLSCStrategy::getRequestedParameters()
const
520 Teuchos::RCP<Teuchos::ParameterList> result;
521 Teuchos::RCP<Teuchos::ParameterList> pl = rcp(
new Teuchos::ParameterList());
524 RCP<Teuchos::ParameterList> fList = invFactoryF_->getRequestedParameters();
525 if(fList!=Teuchos::null) {
526 Teuchos::ParameterList::ConstIterator itr;
527 for(itr=fList->begin();itr!=fList->end();++itr)
528 pl->setEntry(itr->first,itr->second);
533 RCP<Teuchos::ParameterList> sList = invFactoryS_->getRequestedParameters();
534 if(sList!=Teuchos::null) {
535 Teuchos::ParameterList::ConstIterator itr;
536 for(itr=sList->begin();itr!=sList->end();++itr)
537 pl->setEntry(itr->first,itr->second);
543 pl->set<Teko::LinearOp>(
"W-Scaling Vector", Teuchos::null,
"W-Scaling Vector");
551bool InvLSCStrategy::updateRequestedParameters(
const Teuchos::ParameterList & pl)
553 Teko_DEBUG_SCOPE(
"InvLSCStrategy::updateRequestedParameters",10);
557 result &= invFactoryF_->updateRequestedParameters(pl);
558 result &= invFactoryS_->updateRequestedParameters(pl);
562 Teko::MultiVector wScale = pl.get<Teko::MultiVector>(
"W-Scaling Vector");
564 if(wScale==Teuchos::null)
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.
@ NotDiag
For user convenience, if Teko recieves this value, exceptions will be thrown.
@ Diagonal
Specifies that just the diagonal is used.
MultiVector getBlock(int i, const BlockedMultiVector &bmv)
Get the ith block from a BlockedMultiVector object.
LinearOp adjoint(ModifiableLinearOp &a)
Construct an implicit adjoint of the linear operators.
An implementation of a state object for block preconditioners.
Preconditioner state for the LSC factory.
LinearOp invMass_
Inverse mass operator ( )
virtual bool isInitialized() const
virtual void addInverse(const std::string &name, const Teko::InverseLinearOp &ilo)
Add a named inverse to the state object.
virtual void setInitialized(bool init=true)
virtual Teko::InverseLinearOp getInverse(const std::string &name) const
Get a named inverse from the state object.