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"
79 using Teuchos::rcp_const_cast;
80 using Teuchos::rcp_dynamic_cast;
91 InvLSCStrategy::InvLSCStrategy()
92 : massMatrix_(Teuchos::null),
93 invFactoryF_(Teuchos::null),
94 invFactoryS_(Teuchos::null),
96 rowZeroingNeeded_(false),
101 scaleType_(Diagonal),
103 assumeStable_(false) {}
105 InvLSCStrategy::InvLSCStrategy(
const Teuchos::RCP<InverseFactory>& factory,
bool rzn)
106 : massMatrix_(Teuchos::null),
107 invFactoryF_(factory),
108 invFactoryS_(factory),
110 rowZeroingNeeded_(rzn),
115 scaleType_(Diagonal),
117 assumeStable_(false) {}
119 InvLSCStrategy::InvLSCStrategy(
const Teuchos::RCP<InverseFactory>& invFactF,
120 const Teuchos::RCP<InverseFactory>& invFactS,
bool rzn)
121 : massMatrix_(Teuchos::null),
122 invFactoryF_(invFactF),
123 invFactoryS_(invFactS),
125 rowZeroingNeeded_(rzn),
130 scaleType_(Diagonal),
132 assumeStable_(false) {}
134 InvLSCStrategy::InvLSCStrategy(
const Teuchos::RCP<InverseFactory>& factory, LinearOp& mass,
137 invFactoryF_(factory),
138 invFactoryS_(factory),
140 rowZeroingNeeded_(rzn),
145 scaleType_(Diagonal),
147 assumeStable_(false) {}
149 InvLSCStrategy::InvLSCStrategy(
const Teuchos::RCP<InverseFactory>& invFactF,
150 const Teuchos::RCP<InverseFactory>& invFactS, LinearOp& mass,
153 invFactoryF_(invFactF),
154 invFactoryS_(invFactS),
156 rowZeroingNeeded_(rzn),
161 scaleType_(Diagonal),
163 assumeStable_(false) {}
168 Teko_DEBUG_SCOPE(
"InvLSCStrategy::buildState", 10);
171 TEUCHOS_ASSERT(lscState != 0);
175 Teko_DEBUG_EXPR(Teuchos::Time timer(
""));
179 Teko_DEBUG_SCOPE(
"LSC::buildState constructing operators", 1);
180 Teko_DEBUG_EXPR(timer.start(
true));
182 initializeState(A, lscState);
184 Teko_DEBUG_EXPR(timer.stop());
185 Teko_DEBUG_MSG(
"LSC::buildState BuildOpsTime = " << timer.totalElapsedTime(), 1);
190 Teko_DEBUG_SCOPE(
"LSC::buildState calculating inverses", 1);
191 Teko_DEBUG_EXPR(timer.start(
true));
193 computeInverses(A, lscState);
195 Teko_DEBUG_EXPR(timer.stop());
196 Teko_DEBUG_MSG(
"LSC::buildState BuildInvTime = " << timer.totalElapsedTime(), 1);
202 LinearOp InvLSCStrategy::getInvBQBt(
const BlockedLinearOp& ,
207 LinearOp InvLSCStrategy::getInvBHBt(
const BlockedLinearOp& ,
212 LinearOp InvLSCStrategy::getInvF(
const BlockedLinearOp& ,
217 LinearOp InvLSCStrategy::getOuterStabilization(
const BlockedLinearOp& ,
220 TEUCHOS_ASSERT(lscState != 0);
223 return lscState->
aiD_;
226 LinearOp InvLSCStrategy::getInvMass(
const BlockedLinearOp& ,
229 TEUCHOS_ASSERT(lscState != 0);
235 LinearOp InvLSCStrategy::getHScaling(
const BlockedLinearOp& A,
237 if (hScaling_ != Teuchos::null)
return hScaling_;
238 return getInvMass(A, state);
242 void InvLSCStrategy::initializeState(
const BlockedLinearOp& A,
LSCPrecondState* state)
const {
243 Teko_DEBUG_SCOPE(
"InvLSCStrategy::initializeState", 10);
245 const LinearOp F = getBlock(0, 0, A);
246 const LinearOp Bt = getBlock(0, 1, A);
247 const LinearOp B = getBlock(1, 0, A);
248 const LinearOp C = getBlock(1, 1, A);
251 LinearOp G = isSymmetric_ ? Bt : adjoint(D);
253 bool isStabilized = assumeStable_ ?
false : (not isZeroOp(C));
261 if (massMatrix_ == Teuchos::null) {
263 "LSC::initializeState Build Scaling <F> type \"" << getDiagonalName(scaleType_) <<
"\"", 1);
264 state->
invMass_ = getInvDiagonalOp(F, scaleType_);
265 }
else if (state->
invMass_ == Teuchos::null) {
267 "LSC::initializeState Build Scaling <mass> type \"" << getDiagonalName(scaleType_) <<
"\"",
269 state->
invMass_ = getInvDiagonalOp(massMatrix_, scaleType_);
275 Teko_DEBUG_MSG(
"Computed BQBt", 10);
278 if (wScaling_ != Teuchos::null && hScaling_ == Teuchos::null) {
280 RCP<const Thyra::VectorBase<double> > w = wScaling_->col(0);
281 RCP<const Thyra::VectorBase<double> > iQu =
282 rcp_dynamic_cast<
const Thyra::DiagonalLinearOpBase<double> >(state->
invMass_)->getDiag();
283 RCP<Thyra::VectorBase<double> > h = Thyra::createMember(iQu->space());
285 Thyra::put_scalar(0.0, h.ptr());
286 Thyra::ele_wise_prod(1.0, *w, *iQu, h.ptr());
287 hScaling_ = Teuchos::rcp(
new Thyra::DefaultDiagonalLinearOp<double>(h));
290 LinearOp H = hScaling_;
291 if (H == Teuchos::null && not isSymmetric_) H = state->
invMass_;
294 if (H == Teuchos::null)
297 RCP<Teuchos::Time> time =
298 Teuchos::TimeMonitor::getNewTimer(
"InvLSCStrategy::initializeState Build BHBt");
299 Teuchos::TimeMonitor timer(*time);
302 state->
BHBt_ = explicitMultiply(D, H, G, state->
BHBt_);
306 if (not isStabilized) {
311 state->
aiD_ = Teuchos::null;
320 if (!Teko::TpetraHelpers::isTpetraLinearOp(F)) {
321 #ifdef TEKO_HAVE_EPETRA
322 const RCP<const Epetra_Operator> epF = Thyra::get_Epetra_Operator(*F);
323 if (epF != Teuchos::null && rowZeroingNeeded_) {
325 const RCP<const Epetra_CrsMatrix> crsF = rcp_dynamic_cast<
const Epetra_CrsMatrix>(epF);
328 if (crsF != Teuchos::null) {
329 std::vector<int> zeroIndices;
332 Teko::Epetra::identityRowIndices(crsF->RowMap(), *crsF, zeroIndices);
339 throw std::logic_error(
340 "InvLSCStrategy::initializeState is trying to use "
341 "Epetra code, but TEKO is not built with Epetra!");
346 RCP<const Tpetra::CrsMatrix<ST, LO, GO, NT> > crsF =
347 Teko::TpetraHelpers::getTpetraCrsMatrix(F, &scalar, &transp);
349 std::vector<GO> zeroIndices;
352 Teko::TpetraHelpers::identityRowIndices(*crsF->getRowMap(), *crsF, zeroIndices);
355 modF = Thyra::tpetraLinearOp<ST, LO, GO, NT>(
356 Thyra::tpetraVectorSpace<ST, LO, GO, NT>(crsF->getDomainMap()),
357 Thyra::tpetraVectorSpace<ST, LO, GO, NT>(crsF->getRangeMap()),
362 Teko_DEBUG_MSG(
"Calculating gamma", 10);
363 LinearOp iQuF = multiply(state->
invMass_, modF);
366 Teko::LinearOp stabMatrix;
367 state->
gamma_ = std::fabs(Teko::computeSpectralRad(iQuF, 5e-2,
false, eigSolveParam_)) / 3.0;
368 Teko_DEBUG_MSG(
"Calculated gamma", 10);
369 if (userPresStabMat_ != Teuchos::null) {
370 Teko::LinearOp invDGl = Teko::getInvDiagonalOp(userPresStabMat_);
371 Teko::LinearOp gammaOp = multiply(invDGl, C);
372 state->
gamma_ *= std::fabs(Teko::computeSpectralRad(gammaOp, 5e-2,
false, eigSolveParam_));
373 stabMatrix = userPresStabMat_;
379 LinearOp invDiagF = getInvDiagonalOp(F);
380 Teko::ModifiableLinearOp modB_idF_Bt = state->
getInverse(
"BidFBt");
381 modB_idF_Bt = explicitMultiply(B, invDiagF, Bt, modB_idF_Bt);
383 const LinearOp B_idF_Bt = modB_idF_Bt;
385 MultiVector vec_D = getDiagonal(B_idF_Bt);
386 update(-1.0, getDiagonal(C), 1.0, vec_D);
387 const LinearOp invD = buildInvDiagonal(vec_D,
"inv(D)");
389 Teko_DEBUG_MSG(
"Calculating alpha", 10);
390 const LinearOp BidFBtidD = multiply<double>(B_idF_Bt, invD);
391 double num = std::fabs(Teko::computeSpectralRad(BidFBtidD, 5e-2,
false, eigSolveParam_));
392 Teko_DEBUG_MSG(
"Calculated alpha", 10);
393 state->
alpha_ = 1.0 / num;
394 state->
aiD_ = Thyra::scale(state->
alpha_, invD);
397 Teko::ModifiableLinearOp BQBtmC = state->
getInverse(
"BQBtmC");
398 BQBtmC = explicitAdd(state->
BQBt_, scale(-state->
gamma_, stabMatrix), BQBtmC);
402 Teko::ModifiableLinearOp BHBtmC = state->
getInverse(
"BHBtmC");
403 if (H == Teuchos::null)
406 BHBtmC = explicitAdd(state->
BHBt_, scale(-state->
gamma_, stabMatrix), BHBtmC);
410 Teko_DEBUG_MSG_BEGIN(5) DEBUG_STREAM <<
"LSC Gamma Parameter = " << state->
gamma_ << std::endl;
411 DEBUG_STREAM <<
"LSC Alpha Parameter = " << state->
alpha_ << std::endl;
422 void InvLSCStrategy::computeInverses(
const BlockedLinearOp& A,
LSCPrecondState* state)
const {
423 Teko_DEBUG_SCOPE(
"InvLSCStrategy::computeInverses", 10);
424 Teko_DEBUG_EXPR(Teuchos::Time invTimer(
""));
426 const LinearOp F = getBlock(0, 0, A);
431 Teko_DEBUG_MSG(
"LSC::computeInverses Building inv(F)", 1);
432 Teko_DEBUG_EXPR(invTimer.start(
true));
433 InverseLinearOp invF = state->
getInverse(
"invF");
434 if (invF == Teuchos::null) {
440 Teko_DEBUG_EXPR(invTimer.stop());
441 Teko_DEBUG_MSG(
"LSC::computeInverses GetInvF = " << invTimer.totalElapsedTime(), 1);
446 Teko_DEBUG_MSG(
"LSC::computeInverses Building inv(BQBtmC)", 1);
447 Teko_DEBUG_EXPR(invTimer.start(
true));
448 const LinearOp BQBt = state->
getInverse(
"BQBtmC");
449 InverseLinearOp invBQBt = state->
getInverse(
"invBQBtmC");
450 if (invBQBt == Teuchos::null) {
456 Teko_DEBUG_EXPR(invTimer.stop());
457 Teko_DEBUG_MSG(
"LSC::computeInverses GetInvBQBt = " << invTimer.totalElapsedTime(), 1);
462 ModifiableLinearOp invBHBt = state->
getInverse(
"invBHBtmC");
463 if (hScaling_ != Teuchos::null || not isSymmetric_) {
465 Teko_DEBUG_MSG(
"LSC::computeInverses Building inv(BHBtmC)", 1);
466 Teko_DEBUG_EXPR(invTimer.start(
true));
467 const LinearOp BHBt = state->
getInverse(
"BHBtmC");
468 if (invBHBt == Teuchos::null) {
474 Teko_DEBUG_EXPR(invTimer.stop());
475 Teko_DEBUG_MSG(
"LSC::computeInverses GetInvBHBt = " << invTimer.totalElapsedTime(), 1);
476 }
else if (invBHBt == Teuchos::null) {
483 void InvLSCStrategy::initializeFromParameterList(
const Teuchos::ParameterList& pl,
484 const InverseLibrary& invLib) {
486 std::string invStr =
"", invVStr =
"", invPStr =
"";
487 bool rowZeroing =
true;
489 scaleType_ = Diagonal;
492 if (pl.isParameter(
"Inverse Type")) invStr = pl.get<std::string>(
"Inverse Type");
493 if (pl.isParameter(
"Inverse Velocity Type"))
494 invVStr = pl.get<std::string>(
"Inverse Velocity Type");
495 if (pl.isParameter(
"Inverse Pressure Type"))
496 invPStr = pl.get<std::string>(
"Inverse Pressure Type");
497 if (pl.isParameter(
"Ignore Boundary Rows")) rowZeroing = pl.get<
bool>(
"Ignore Boundary Rows");
498 if (pl.isParameter(
"Use LDU")) useLDU = pl.get<
bool>(
"Use LDU");
499 if (pl.isParameter(
"Use Mass Scaling")) useMass_ = pl.get<
bool>(
"Use Mass Scaling");
502 if (pl.isParameter(
"Use W-Scaling")) useWScaling_ = pl.get<
bool>(
"Use W-Scaling");
503 if (pl.isParameter(
"Eigen Solver Iterations"))
504 eigSolveParam_ = pl.get<
int>(
"Eigen Solver Iterations");
505 if (pl.isParameter(
"Scaling Type")) {
506 scaleType_ = getDiagonalType(pl.get<std::string>(
"Scaling Type"));
507 TEUCHOS_TEST_FOR_EXCEPT(scaleType_ == NotDiag);
509 if (pl.isParameter(
"Assume Stable Discretization"))
510 assumeStable_ = pl.get<
bool>(
"Assume Stable Discretization");
512 Teko_DEBUG_MSG_BEGIN(5) DEBUG_STREAM <<
"LSC Inverse Strategy Parameters: " << std::endl;
513 DEBUG_STREAM <<
" inv type = \"" << invStr <<
"\"" << std::endl;
514 DEBUG_STREAM <<
" inv v type = \"" << invVStr <<
"\"" << std::endl;
515 DEBUG_STREAM <<
" inv p type = \"" << invPStr <<
"\"" << std::endl;
516 DEBUG_STREAM <<
" bndry rows = " << rowZeroing << std::endl;
517 DEBUG_STREAM <<
" use ldu = " << useLDU << std::endl;
518 DEBUG_STREAM <<
" use mass = " << useMass_ << std::endl;
519 DEBUG_STREAM <<
" use w-scaling = " << useWScaling_ << std::endl;
520 DEBUG_STREAM <<
" assume stable = " << assumeStable_ << std::endl;
521 DEBUG_STREAM <<
" scale type = " << getDiagonalName(scaleType_) << std::endl;
522 DEBUG_STREAM <<
"LSC Inverse Strategy Parameter list: " << std::endl;
523 pl.print(DEBUG_STREAM);
527 #if defined(Teko_ENABLE_Amesos)
528 if (invStr ==
"") invStr =
"Amesos";
529 #elif defined(Teko_ENABLE_Amesos2)
530 if (invStr ==
"") invStr =
"Amesos2";
532 if (invVStr ==
"") invVStr = invStr;
533 if (invPStr ==
"") invPStr = invStr;
536 invFactoryF_ = invLib.getInverseFactory(invVStr);
537 invFactoryS_ = invFactoryF_;
538 if (invVStr != invPStr)
539 invFactoryS_ = invLib.getInverseFactory(invPStr);
542 setUseFullLDU(useLDU);
543 setRowZeroing(rowZeroing);
546 Teuchos::RCP<Teko::RequestHandler> rh = getRequestHandler();
547 rh->preRequest<Teko::LinearOp>(Teko::RequestMesg(
"Velocity Mass Matrix"));
548 Teko::LinearOp mass = rh->request<Teko::LinearOp>(Teko::RequestMesg(
"Velocity Mass Matrix"));
554 Teuchos::RCP<Teuchos::ParameterList> InvLSCStrategy::getRequestedParameters()
const {
555 Teuchos::RCP<Teuchos::ParameterList> result;
556 Teuchos::RCP<Teuchos::ParameterList> pl = rcp(
new Teuchos::ParameterList());
559 RCP<Teuchos::ParameterList> fList = invFactoryF_->getRequestedParameters();
560 if (fList != Teuchos::null) {
561 Teuchos::ParameterList::ConstIterator itr;
562 for (itr = fList->begin(); itr != fList->end(); ++itr) pl->setEntry(itr->first, itr->second);
567 RCP<Teuchos::ParameterList> sList = invFactoryS_->getRequestedParameters();
568 if (sList != Teuchos::null) {
569 Teuchos::ParameterList::ConstIterator itr;
570 for (itr = sList->begin(); itr != sList->end(); ++itr) pl->setEntry(itr->first, itr->second);
576 pl->set<Teko::LinearOp>(
"W-Scaling Vector", Teuchos::null,
"W-Scaling Vector");
584 bool InvLSCStrategy::updateRequestedParameters(
const Teuchos::ParameterList& pl) {
585 Teko_DEBUG_SCOPE(
"InvLSCStrategy::updateRequestedParameters", 10);
589 result &= invFactoryF_->updateRequestedParameters(pl);
590 result &= invFactoryS_->updateRequestedParameters(pl);
594 Teko::MultiVector wScale = pl.get<Teko::MultiVector>(
"W-Scaling Vector");
596 if (wScale == Teuchos::null)
void rebuildInverse(const InverseFactory &factory, const LinearOp &A, InverseLinearOp &invA)
virtual Teko::InverseLinearOp getInverse(const std::string &name) const
Get a named inverse from the state object.
An implementation of a state object for block preconditioners.
virtual bool isInitialized() const
virtual void addInverse(const std::string &name, const Teko::InverseLinearOp &ilo)
Add a named inverse to the state object.
LinearOp invMass_
Inverse mass operator ( )
Preconditioner state for the LSC factory.
InverseLinearOp buildInverse(const InverseFactory &factory, const LinearOp &A)
Build an inverse operator using a factory and a linear operator.
virtual void setInitialized(bool init=true)