Teko Version of the Day
Teko_InvLSCStrategy.cpp
1/*
2// @HEADER
3//
4// ***********************************************************************
5//
6// Teko: A package for block and physics based preconditioning
7// Copyright 2010 Sandia Corporation
8//
9// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
10// the U.S. Government retains certain rights in this software.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric C. Cyr (eccyr@sandia.gov)
40//
41// ***********************************************************************
42//
43// @HEADER
44
45*/
46
47#include "NS/Teko_InvLSCStrategy.hpp"
48
49#include "Thyra_DefaultDiagonalLinearOp.hpp"
50#include "Thyra_EpetraThyraWrappers.hpp"
51#include "Thyra_get_Epetra_Operator.hpp"
52#include "Thyra_EpetraLinearOp.hpp"
53#include "Thyra_VectorStdOps.hpp"
54
55#include "Epetra_Vector.h"
56#include "Epetra_Map.h"
57
58#include "EpetraExt_RowMatrixOut.h"
59#include "EpetraExt_MultiVectorOut.h"
60#include "EpetraExt_VectorOut.h"
61
62#include "Teuchos_Time.hpp"
63#include "Teuchos_TimeMonitor.hpp"
64
65// Teko includes
66#include "Teko_Utilities.hpp"
67#include "NS/Teko_LSCPreconditionerFactory.hpp"
68#include "Teko_EpetraHelpers.hpp"
69#include "Teko_EpetraOperatorWrapper.hpp"
70#include "Teko_TpetraHelpers.hpp"
71
72#include "Thyra_TpetraLinearOp.hpp"
73
74using Teuchos::RCP;
75using Teuchos::rcp_dynamic_cast;
76using Teuchos::rcp_const_cast;
77
78namespace Teko {
79namespace NS {
80
82// InvLSCStrategy Implementation
84
85// constructors
87InvLSCStrategy::InvLSCStrategy()
88 : massMatrix_(Teuchos::null), invFactoryF_(Teuchos::null), invFactoryS_(Teuchos::null), eigSolveParam_(5)
89 , rowZeroingNeeded_(false), useFullLDU_(false), useMass_(false), useLumping_(false), useWScaling_(false), scaleType_(Diagonal)
90 , isSymmetric_(true), assumeStable_(false)
91{ }
92
93InvLSCStrategy::InvLSCStrategy(const Teuchos::RCP<InverseFactory> & factory,bool rzn)
94 : massMatrix_(Teuchos::null), invFactoryF_(factory), invFactoryS_(factory), eigSolveParam_(5), rowZeroingNeeded_(rzn)
95 , useFullLDU_(false), useMass_(false), useLumping_(false), useWScaling_(false), scaleType_(Diagonal)
96 , isSymmetric_(true), assumeStable_(false)
97{ }
98
99InvLSCStrategy::InvLSCStrategy(const Teuchos::RCP<InverseFactory> & invFactF,
100 const Teuchos::RCP<InverseFactory> & invFactS,
101 bool rzn)
102 : massMatrix_(Teuchos::null), invFactoryF_(invFactF), invFactoryS_(invFactS), eigSolveParam_(5), rowZeroingNeeded_(rzn)
103 , useFullLDU_(false), useMass_(false), useLumping_(false), useWScaling_(false), scaleType_(Diagonal)
104 , isSymmetric_(true), assumeStable_(false)
105{ }
106
107InvLSCStrategy::InvLSCStrategy(const Teuchos::RCP<InverseFactory> & factory,LinearOp & mass,bool rzn)
108 : massMatrix_(mass), invFactoryF_(factory), invFactoryS_(factory), eigSolveParam_(5), rowZeroingNeeded_(rzn)
109 , useFullLDU_(false), useMass_(false), useLumping_(false), useWScaling_(false), scaleType_(Diagonal)
110 , isSymmetric_(true), assumeStable_(false)
111{ }
112
113InvLSCStrategy::InvLSCStrategy(const Teuchos::RCP<InverseFactory> & invFactF,
114 const Teuchos::RCP<InverseFactory> & invFactS,
115 LinearOp & mass,bool rzn)
116 : massMatrix_(mass), invFactoryF_(invFactF), invFactoryS_(invFactS), eigSolveParam_(5), rowZeroingNeeded_(rzn)
117 , useFullLDU_(false), useMass_(false), useLumping_(false), useWScaling_(false), scaleType_(Diagonal)
118 , isSymmetric_(true), assumeStable_(false)
119{ }
120
122
123void InvLSCStrategy::buildState(BlockedLinearOp & A,BlockPreconditionerState & state) const
124{
125 Teko_DEBUG_SCOPE("InvLSCStrategy::buildState",10);
126
127 LSCPrecondState * lscState = dynamic_cast<LSCPrecondState*>(&state);
128 TEUCHOS_ASSERT(lscState!=0);
129
130 // if neccessary save state information
131 if(not lscState->isInitialized()) {
132 Teko_DEBUG_EXPR(Teuchos::Time timer(""));
133
134 // construct operators
135 {
136 Teko_DEBUG_SCOPE("LSC::buildState constructing operators",1);
137 Teko_DEBUG_EXPR(timer.start(true));
138
139 initializeState(A,lscState);
140
141 Teko_DEBUG_EXPR(timer.stop());
142 Teko_DEBUG_MSG("LSC::buildState BuildOpsTime = " << timer.totalElapsedTime(),1);
143 }
144
145 // Build the inverses
146 {
147 Teko_DEBUG_SCOPE("LSC::buildState calculating inverses",1);
148 Teko_DEBUG_EXPR(timer.start(true));
149
150 computeInverses(A,lscState);
151
152 Teko_DEBUG_EXPR(timer.stop());
153 Teko_DEBUG_MSG("LSC::buildState BuildInvTime = " << timer.totalElapsedTime(),1);
154 }
155 }
156}
157
158// functions inherited from LSCStrategy
159LinearOp InvLSCStrategy::getInvBQBt(const BlockedLinearOp & /* A */,BlockPreconditionerState & state) const
160{
161 return state.getInverse("invBQBtmC");
162}
163
164LinearOp InvLSCStrategy::getInvBHBt(const BlockedLinearOp & /* A */,BlockPreconditionerState & state) const
165{
166 return state.getInverse("invBHBtmC");
167}
168
169LinearOp InvLSCStrategy::getInvF(const BlockedLinearOp & /* A */,BlockPreconditionerState & state) const
170{
171 return state.getInverse("invF");
172}
173
174LinearOp InvLSCStrategy::getOuterStabilization(const BlockedLinearOp & /* A */,BlockPreconditionerState & state) const
175{
176 LSCPrecondState * lscState = dynamic_cast<LSCPrecondState*>(&state);
177 TEUCHOS_ASSERT(lscState!=0);
178 TEUCHOS_ASSERT(lscState->isInitialized())
179
180 return lscState->aiD_;
181}
182
183LinearOp InvLSCStrategy::getInvMass(const BlockedLinearOp & /* A */,BlockPreconditionerState & state) const
184{
185 LSCPrecondState * lscState = dynamic_cast<LSCPrecondState*>(&state);
186 TEUCHOS_ASSERT(lscState!=0);
187 TEUCHOS_ASSERT(lscState->isInitialized())
188
189 return lscState->invMass_;
190}
191
192LinearOp InvLSCStrategy::getHScaling(const BlockedLinearOp & A,BlockPreconditionerState & state) const
193{
194 if(hScaling_!=Teuchos::null) return hScaling_;
195 return getInvMass(A,state);
196}
197
199void InvLSCStrategy::initializeState(const BlockedLinearOp & A,LSCPrecondState * state) const
200{
201 Teko_DEBUG_SCOPE("InvLSCStrategy::initializeState",10);
202
203 const LinearOp F = getBlock(0,0,A);
204 const LinearOp Bt = getBlock(0,1,A);
205 const LinearOp B = getBlock(1,0,A);
206 const LinearOp C = getBlock(1,1,A);
207
208 LinearOp D = B;
209 LinearOp G = isSymmetric_ ? Bt : adjoint(D);
210
211 bool isStabilized = assumeStable_ ? false : (not isZeroOp(C));
212
213 // The logic follows like this
214 // if there is no mass matrix available --> build from F
215 // if there is a mass matrix and the inverse hasn't yet been built
216 // --> build from the mass matrix
217 // otherwise, there is already an invMass_ matrix that is appropriate
218 // --> use that one
219 if(massMatrix_==Teuchos::null) {
220 Teko_DEBUG_MSG("LSC::initializeState Build Scaling <F> type \""
221 << getDiagonalName(scaleType_) << "\"" ,1);
222 state->invMass_ = getInvDiagonalOp(F,scaleType_);
223 }
224 else if(state->invMass_==Teuchos::null) {
225 Teko_DEBUG_MSG("LSC::initializeState Build Scaling <mass> type \""
226 << getDiagonalName(scaleType_) << "\"" ,1);
227 state->invMass_ = getInvDiagonalOp(massMatrix_,scaleType_);
228 }
229 // else "invMass_" should be set and there is no reason to rebuild it
230
231 // compute BQBt
232 state->BQBt_ = explicitMultiply(B,state->invMass_,Bt,state->BQBt_);
233 Teko_DEBUG_MSG("Computed BQBt",10);
234
235 // if there is no H-Scaling
236 if(wScaling_!=Teuchos::null && hScaling_==Teuchos::null) {
237 // from W vector build H operator scaling
238 RCP<const Thyra::VectorBase<double> > w = wScaling_->col(0);
239 RCP<const Thyra::VectorBase<double> > iQu
240 = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<double> >(state->invMass_)->getDiag();
241 RCP<Thyra::VectorBase<double> > h = Thyra::createMember(iQu->space());
242
243 Thyra::put_scalar(0.0,h.ptr());
244 Thyra::ele_wise_prod(1.0,*w,*iQu,h.ptr());
245 hScaling_ = Teuchos::rcp(new Thyra::DefaultDiagonalLinearOp<double>(h));
246 }
247
248 LinearOp H = hScaling_;
249 if(H==Teuchos::null && not isSymmetric_)
250 H = state->invMass_;
251
252 // setup the scaling operator
253 if(H==Teuchos::null)
254 state->BHBt_ = state->BQBt_;
255 else {
256 RCP<Teuchos::Time> time = Teuchos::TimeMonitor::getNewTimer("InvLSCStrategy::initializeState Build BHBt");
257 Teuchos::TimeMonitor timer(*time);
258
259 // compute BHBt
260 state->BHBt_ = explicitMultiply(D,H,G,state->BHBt_);
261 }
262
263 // if this is a stable discretization...we are done!
264 if(not isStabilized) {
265 state->addInverse("BQBtmC",state->BQBt_);
266 state->addInverse("BHBtmC",state->BHBt_);
267 state->gamma_ = 0.0;
268 state->alpha_ = 0.0;
269 state->aiD_ = Teuchos::null;
270
271 state->setInitialized(true);
272
273 return;
274 }
275
276 // for Epetra_CrsMatrix...zero out certain rows: this ensures spectral radius is correct
277 LinearOp modF = F;
278 if(!Teko::TpetraHelpers::isTpetraLinearOp(F)){ // Epetra
279 const RCP<const Epetra_Operator> epF = Thyra::get_Epetra_Operator(*F);
280 if(epF!=Teuchos::null && rowZeroingNeeded_) {
281 // try to get a CRS matrix
282 const RCP<const Epetra_CrsMatrix> crsF = rcp_dynamic_cast<const Epetra_CrsMatrix>(epF);
283
284 // if it is a CRS matrix get rows that need to be zeroed
285 if(crsF!=Teuchos::null) {
286 std::vector<int> zeroIndices;
287
288 // get rows in need of zeroing
289 Teko::Epetra::identityRowIndices(crsF->RowMap(), *crsF,zeroIndices);
290
291 // build an operator that zeros those rows
292 modF = Thyra::epetraLinearOp(rcp(new Teko::Epetra::ZeroedOperator(zeroIndices,crsF)));
293 }
294 }
295 } else { //Tpetra
296 ST scalar = 0.0;
297 bool transp = false;
298 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > crsF = Teko::TpetraHelpers::getTpetraCrsMatrix(F, &scalar, &transp);
299
300 std::vector<GO> zeroIndices;
301
302 // get rows in need of zeroing
303 Teko::TpetraHelpers::identityRowIndices(*crsF->getRowMap(), *crsF,zeroIndices);
304
305 // build an operator that zeros those rows
306 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)));
307 }
308
309 // compute gamma
310 Teko_DEBUG_MSG("Calculating gamma",10);
311 LinearOp iQuF = multiply(state->invMass_,modF);
312
313 // do 6 power iterations to compute spectral radius: EHSST2007 Eq. 4.28
314 Teko::LinearOp stabMatrix; // this is the pressure stabilization matrix to use
315 state->gamma_ = std::fabs(Teko::computeSpectralRad(iQuF,5e-2,false,eigSolveParam_))/3.0;
316 Teko_DEBUG_MSG("Calculated gamma",10);
317 if(userPresStabMat_!=Teuchos::null) {
318 Teko::LinearOp invDGl = Teko::getInvDiagonalOp(userPresStabMat_);
319 Teko::LinearOp gammaOp = multiply(invDGl,C);
320 state->gamma_ *= std::fabs(Teko::computeSpectralRad(gammaOp,5e-2,false,eigSolveParam_));
321 stabMatrix = userPresStabMat_;
322 } else
323 stabMatrix = C;
324
325 // compute alpha scaled inv(D): EHSST2007 Eq. 4.29
326 // construct B_idF_Bt and save it for refilling later: This could reuse BQBt graph
327 LinearOp invDiagF = getInvDiagonalOp(F);
328 Teko::ModifiableLinearOp modB_idF_Bt = state->getInverse("BidFBt");
329 modB_idF_Bt = explicitMultiply(B,invDiagF,Bt,modB_idF_Bt);
330 state->addInverse("BidFBt",modB_idF_Bt);
331 const LinearOp B_idF_Bt = modB_idF_Bt;
332
333 MultiVector vec_D = getDiagonal(B_idF_Bt); // this memory could be reused
334 update(-1.0,getDiagonal(C),1.0,vec_D); // vec_D = diag(B*inv(diag(F))*Bt)-diag(C)
335 const LinearOp invD = buildInvDiagonal(vec_D,"inv(D)");
336
337 Teko_DEBUG_MSG("Calculating alpha",10);
338 const LinearOp BidFBtidD = multiply<double>(B_idF_Bt,invD);
339 double num = std::fabs(Teko::computeSpectralRad(BidFBtidD,5e-2,false,eigSolveParam_));
340 Teko_DEBUG_MSG("Calculated alpha",10);
341 state->alpha_ = 1.0/num;
342 state->aiD_ = Thyra::scale(state->alpha_,invD);
343
344 // now build B*Q*Bt-gamma*C
345 Teko::ModifiableLinearOp BQBtmC = state->getInverse("BQBtmC");
346 BQBtmC = explicitAdd(state->BQBt_,scale(-state->gamma_,stabMatrix),BQBtmC);
347 state->addInverse("BQBtmC",BQBtmC);
348
349 // now build B*H*Bt-gamma*C
350 Teko::ModifiableLinearOp BHBtmC = state->getInverse("BHBtmC");
351 if(H==Teuchos::null)
352 BHBtmC = BQBtmC;
353 else {
354 BHBtmC = explicitAdd(state->BHBt_,scale(-state->gamma_,stabMatrix),BHBtmC);
355 }
356 state->addInverse("BHBtmC",BHBtmC);
357
358 Teko_DEBUG_MSG_BEGIN(5)
359 DEBUG_STREAM << "LSC Gamma Parameter = " << state->gamma_ << std::endl;
360 DEBUG_STREAM << "LSC Alpha Parameter = " << state->alpha_ << std::endl;
361 Teko_DEBUG_MSG_END()
362
363 state->setInitialized(true);
364}
365
371void InvLSCStrategy::computeInverses(const BlockedLinearOp & A,LSCPrecondState * state) const
372{
373 Teko_DEBUG_SCOPE("InvLSCStrategy::computeInverses",10);
374 Teko_DEBUG_EXPR(Teuchos::Time invTimer(""));
375
376 const LinearOp F = getBlock(0,0,A);
377
379
380 // (re)build the inverse of F
381 Teko_DEBUG_MSG("LSC::computeInverses Building inv(F)",1);
382 Teko_DEBUG_EXPR(invTimer.start(true));
383 InverseLinearOp invF = state->getInverse("invF");
384 if(invF==Teuchos::null) {
385 invF = buildInverse(*invFactoryF_,F);
386 state->addInverse("invF",invF);
387 } else {
388 rebuildInverse(*invFactoryF_,F,invF);
389 }
390 Teko_DEBUG_EXPR(invTimer.stop());
391 Teko_DEBUG_MSG("LSC::computeInverses GetInvF = " << invTimer.totalElapsedTime(),1);
392
394
395 // (re)build the inverse of BQBt
396 Teko_DEBUG_MSG("LSC::computeInverses Building inv(BQBtmC)",1);
397 Teko_DEBUG_EXPR(invTimer.start(true));
398 const LinearOp BQBt = state->getInverse("BQBtmC");
399 InverseLinearOp invBQBt = state->getInverse("invBQBtmC");
400 if(invBQBt==Teuchos::null) {
401 invBQBt = buildInverse(*invFactoryS_,BQBt);
402 state->addInverse("invBQBtmC",invBQBt);
403 } else {
404 rebuildInverse(*invFactoryS_,BQBt,invBQBt);
405 }
406 Teko_DEBUG_EXPR(invTimer.stop());
407 Teko_DEBUG_MSG("LSC::computeInverses GetInvBQBt = " << invTimer.totalElapsedTime(),1);
408
410
411 // Compute the inverse of BHBt or just use BQBt
412 ModifiableLinearOp invBHBt = state->getInverse("invBHBtmC");
413 if(hScaling_!=Teuchos::null || not isSymmetric_) {
414 // (re)build the inverse of BHBt
415 Teko_DEBUG_MSG("LSC::computeInverses Building inv(BHBtmC)",1);
416 Teko_DEBUG_EXPR(invTimer.start(true));
417 const LinearOp BHBt = state->getInverse("BHBtmC");
418 if(invBHBt==Teuchos::null) {
419 invBHBt = buildInverse(*invFactoryS_,BHBt);
420 state->addInverse("invBHBtmC",invBHBt);
421 } else {
422 rebuildInverse(*invFactoryS_,BHBt,invBHBt);
423 }
424 Teko_DEBUG_EXPR(invTimer.stop());
425 Teko_DEBUG_MSG("LSC::computeInverses GetInvBHBt = " << invTimer.totalElapsedTime(),1);
426 }
427 else if(invBHBt==Teuchos::null) {
428 // just use the Q version
429 state->addInverse("invBHBtmC",invBQBt);
430 }
431}
432
434void InvLSCStrategy::initializeFromParameterList(const Teuchos::ParameterList & pl,const InverseLibrary & invLib)
435{
436 // get string specifying inverse
437 std::string invStr="", invVStr="", invPStr="";
438 bool rowZeroing = true;
439 bool useLDU = false;
440 scaleType_ = Diagonal;
441
442 // "parse" the parameter list
443 if(pl.isParameter("Inverse Type"))
444 invStr = pl.get<std::string>("Inverse Type");
445 if(pl.isParameter("Inverse Velocity Type"))
446 invVStr = pl.get<std::string>("Inverse Velocity Type");
447 if(pl.isParameter("Inverse Pressure Type"))
448 invPStr = pl.get<std::string>("Inverse Pressure Type");
449 if(pl.isParameter("Ignore Boundary Rows"))
450 rowZeroing = pl.get<bool>("Ignore Boundary Rows");
451 if(pl.isParameter("Use LDU"))
452 useLDU = pl.get<bool>("Use LDU");
453 if(pl.isParameter("Use Mass Scaling"))
454 useMass_ = pl.get<bool>("Use Mass Scaling");
455 // if(pl.isParameter("Use Lumping"))
456 // useLumping_ = pl.get<bool>("Use Lumping");
457 if(pl.isParameter("Use W-Scaling"))
458 useWScaling_ = pl.get<bool>("Use W-Scaling");
459 if(pl.isParameter("Eigen Solver Iterations"))
460 eigSolveParam_ = pl.get<int>("Eigen Solver Iterations");
461 if(pl.isParameter("Scaling Type")) {
462 scaleType_ = getDiagonalType(pl.get<std::string>("Scaling Type"));
463 TEUCHOS_TEST_FOR_EXCEPT(scaleType_==NotDiag);
464 }
465 if(pl.isParameter("Assume Stable Discretization"))
466 assumeStable_ = pl.get<bool>("Assume Stable Discretization");
467
468 Teko_DEBUG_MSG_BEGIN(5)
469 DEBUG_STREAM << "LSC Inverse Strategy Parameters: " << std::endl;
470 DEBUG_STREAM << " inv type = \"" << invStr << "\"" << std::endl;
471 DEBUG_STREAM << " inv v type = \"" << invVStr << "\"" << std::endl;
472 DEBUG_STREAM << " inv p type = \"" << invPStr << "\"" << std::endl;
473 DEBUG_STREAM << " bndry rows = " << rowZeroing << std::endl;
474 DEBUG_STREAM << " use ldu = " << useLDU << std::endl;
475 DEBUG_STREAM << " use mass = " << useMass_ << std::endl;
476 DEBUG_STREAM << " use w-scaling = " << useWScaling_ << std::endl;
477 DEBUG_STREAM << " assume stable = " << assumeStable_ << std::endl;
478 DEBUG_STREAM << " scale type = " << getDiagonalName(scaleType_) << std::endl;
479 DEBUG_STREAM << "LSC Inverse Strategy Parameter list: " << std::endl;
480 pl.print(DEBUG_STREAM);
481 Teko_DEBUG_MSG_END()
482
483 // set defaults as needed
484 if(invStr=="") invStr = "Amesos";
485 if(invVStr=="") invVStr = invStr;
486 if(invPStr=="") invPStr = invStr;
487
488 // build velocity inverse factory
489 invFactoryF_ = invLib.getInverseFactory(invVStr);
490 invFactoryS_ = invFactoryF_; // by default these are the same
491 if(invVStr!=invPStr) // if different, build pressure inverse factory
492 invFactoryS_ = invLib.getInverseFactory(invPStr);
493
494 // set other parameters
495 setUseFullLDU(useLDU);
496 setRowZeroing(rowZeroing);
497
498 if(useMass_) {
499 Teuchos::RCP<Teko::RequestHandler> rh = getRequestHandler();
500 rh->preRequest<Teko::LinearOp>(Teko::RequestMesg("Velocity Mass Matrix"));
501 Teko::LinearOp mass
502 = rh->request<Teko::LinearOp>(Teko::RequestMesg("Velocity Mass Matrix"));
503 setMassMatrix(mass);
504 }
505
506}
507
509Teuchos::RCP<Teuchos::ParameterList> InvLSCStrategy::getRequestedParameters() const
510{
511 Teuchos::RCP<Teuchos::ParameterList> result;
512 Teuchos::RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList());
513
514 // grab parameters from F solver
515 RCP<Teuchos::ParameterList> fList = invFactoryF_->getRequestedParameters();
516 if(fList!=Teuchos::null) {
517 Teuchos::ParameterList::ConstIterator itr;
518 for(itr=fList->begin();itr!=fList->end();++itr)
519 pl->setEntry(itr->first,itr->second);
520 result = pl;
521 }
522
523 // grab parameters from S solver
524 RCP<Teuchos::ParameterList> sList = invFactoryS_->getRequestedParameters();
525 if(sList!=Teuchos::null) {
526 Teuchos::ParameterList::ConstIterator itr;
527 for(itr=sList->begin();itr!=sList->end();++itr)
528 pl->setEntry(itr->first,itr->second);
529 result = pl;
530 }
531
532 // use the mass matrix
533 if(useWScaling_) {
534 pl->set<Teko::LinearOp>("W-Scaling Vector", Teuchos::null,"W-Scaling Vector");
535 result = pl;
536 }
537
538 return result;
539}
540
542bool InvLSCStrategy::updateRequestedParameters(const Teuchos::ParameterList & pl)
543{
544 Teko_DEBUG_SCOPE("InvLSCStrategy::updateRequestedParameters",10);
545 bool result = true;
546
547 // update requested parameters in solvers
548 result &= invFactoryF_->updateRequestedParameters(pl);
549 result &= invFactoryS_->updateRequestedParameters(pl);
550
551 // use W scaling matrix
552 if(useWScaling_) {
553 Teko::MultiVector wScale = pl.get<Teko::MultiVector>("W-Scaling Vector");
554
555 if(wScale==Teuchos::null)
556 result &= false;
557 else
558 setWScaling(wScale);
559 }
560
561 return result;
562}
563
564} // end namespace NS
565} // end namespace Teko
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.
void rebuildInverse(const InverseFactory &factory, const LinearOp &A, const LinearOp &precOp, InverseLinearOp &invA)
InverseLinearOp buildInverse(const InverseFactory &factory, const LinearOp &A, const LinearOp &precOp)
Preconditioner state for the LSC factory.
LinearOp invMass_
Inverse mass operator ( )
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.