Teko Version of the Day
Teko_InvLSCStrategy.cpp
00001 /*
00002 // @HEADER
00003 // 
00004 // ***********************************************************************
00005 // 
00006 //      Teko: A package for block and physics based preconditioning
00007 //                  Copyright 2010 Sandia Corporation 
00008 //  
00009 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00010 // the U.S. Government retains certain rights in this software.
00011 //  
00012 // Redistribution and use in source and binary forms, with or without
00013 // modification, are permitted provided that the following conditions are
00014 // met:
00015 //  
00016 // 1. Redistributions of source code must retain the above copyright
00017 // notice, this list of conditions and the following disclaimer.
00018 //  
00019 // 2. Redistributions in binary form must reproduce the above copyright
00020 // notice, this list of conditions and the following disclaimer in the
00021 // documentation and/or other materials provided with the distribution.
00022 //  
00023 // 3. Neither the name of the Corporation nor the names of the
00024 // contributors may be used to endorse or promote products derived from
00025 // this software without specific prior written permission. 
00026 //  
00027 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00028 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00029 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00030 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00031 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00032 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00033 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00034 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00035 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 
00036 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00037 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00038 //  
00039 // Questions? Contact Eric C. Cyr (eccyr@sandia.gov)
00040 // 
00041 // ***********************************************************************
00042 // 
00043 // @HEADER
00044 
00045 */
00046 
00047 #include "NS/Teko_InvLSCStrategy.hpp"
00048 
00049 #include "Thyra_DefaultDiagonalLinearOp.hpp"
00050 #include "Thyra_EpetraThyraWrappers.hpp"
00051 #include "Thyra_get_Epetra_Operator.hpp"
00052 #include "Thyra_EpetraLinearOp.hpp"
00053 #include "Thyra_VectorStdOps.hpp"
00054 
00055 #include "Epetra_Vector.h"
00056 #include "Epetra_Map.h"
00057 
00058 #include "EpetraExt_RowMatrixOut.h"
00059 #include "EpetraExt_MultiVectorOut.h"
00060 #include "EpetraExt_VectorOut.h"
00061 
00062 #include "Teuchos_Time.hpp"
00063 #include "Teuchos_TimeMonitor.hpp"
00064 
00065 // Teko includes
00066 #include "Teko_Utilities.hpp"
00067 #include "NS/Teko_LSCPreconditionerFactory.hpp"
00068 #include "Epetra/Teko_EpetraHelpers.hpp"
00069 #include "Epetra/Teko_EpetraOperatorWrapper.hpp"
00070 
00071 using Teuchos::RCP;
00072 using Teuchos::rcp_dynamic_cast;
00073 using Teuchos::rcp_const_cast;
00074 
00075 namespace Teko {
00076 namespace NS {
00077 
00079 // InvLSCStrategy Implementation
00081 
00082 // constructors
00084 InvLSCStrategy::InvLSCStrategy()
00085    : massMatrix_(Teuchos::null), invFactoryF_(Teuchos::null), invFactoryS_(Teuchos::null), eigSolveParam_(5)
00086    , rowZeroingNeeded_(false), useFullLDU_(false), useMass_(false), useLumping_(false), useWScaling_(false), scaleType_(Diagonal)
00087    , isSymmetric_(true)
00088 { }
00089 
00090 InvLSCStrategy::InvLSCStrategy(const Teuchos::RCP<InverseFactory> & factory,bool rzn)
00091    : massMatrix_(Teuchos::null), invFactoryF_(factory), invFactoryS_(factory), eigSolveParam_(5), rowZeroingNeeded_(rzn)
00092    , useFullLDU_(false), useMass_(false), useLumping_(false), useWScaling_(false), scaleType_(Diagonal)
00093    , isSymmetric_(true)
00094 { }
00095 
00096 InvLSCStrategy::InvLSCStrategy(const Teuchos::RCP<InverseFactory> & invFactF,
00097                                const Teuchos::RCP<InverseFactory> & invFactS,
00098                                bool rzn)
00099    : massMatrix_(Teuchos::null), invFactoryF_(invFactF), invFactoryS_(invFactS), eigSolveParam_(5), rowZeroingNeeded_(rzn)
00100    , useFullLDU_(false), useMass_(false), useLumping_(false), useWScaling_(false), scaleType_(Diagonal)
00101    , isSymmetric_(true)
00102 { }
00103 
00104 InvLSCStrategy::InvLSCStrategy(const Teuchos::RCP<InverseFactory> & factory,LinearOp & mass,bool rzn)
00105    : massMatrix_(mass), invFactoryF_(factory), invFactoryS_(factory), eigSolveParam_(5), rowZeroingNeeded_(rzn)
00106    , useFullLDU_(false), useMass_(false), useLumping_(false), useWScaling_(false), scaleType_(Diagonal)
00107    , isSymmetric_(true)
00108 { }
00109 
00110 InvLSCStrategy::InvLSCStrategy(const Teuchos::RCP<InverseFactory> & invFactF,
00111                                const Teuchos::RCP<InverseFactory> & invFactS,
00112                                LinearOp & mass,bool rzn)
00113    : massMatrix_(mass), invFactoryF_(invFactF), invFactoryS_(invFactS), eigSolveParam_(5), rowZeroingNeeded_(rzn)
00114    , useFullLDU_(false), useMass_(false), useLumping_(false), useWScaling_(false), scaleType_(Diagonal)
00115    , isSymmetric_(true)
00116 { }
00117 
00119 
00120 void InvLSCStrategy::buildState(BlockedLinearOp & A,BlockPreconditionerState & state) const
00121 {
00122    Teko_DEBUG_SCOPE("InvLSCStrategy::buildState",10);
00123 
00124    LSCPrecondState * lscState = dynamic_cast<LSCPrecondState*>(&state);
00125    TEUCHOS_ASSERT(lscState!=0);
00126 
00127    // if neccessary save state information
00128    if(not lscState->isInitialized()) {
00129       Teko_DEBUG_EXPR(Teuchos::Time timer(""));
00130 
00131       // construct operators
00132       {
00133          Teko_DEBUG_SCOPE("LSC::buildState constructing operators",1);
00134          Teko_DEBUG_EXPR(timer.start(true));
00135 
00136          initializeState(A,lscState);
00137 
00138          Teko_DEBUG_EXPR(timer.stop());
00139          Teko_DEBUG_MSG("LSC::buildState BuildOpsTime = " << timer.totalElapsedTime(),1);
00140       }
00141 
00142       // Build the inverses
00143       {
00144          Teko_DEBUG_SCOPE("LSC::buildState calculating inverses",1);
00145          Teko_DEBUG_EXPR(timer.start(true));
00146 
00147          computeInverses(A,lscState);
00148 
00149          Teko_DEBUG_EXPR(timer.stop());
00150          Teko_DEBUG_MSG("LSC::buildState BuildInvTime = " << timer.totalElapsedTime(),1);
00151       }
00152    }
00153 }
00154 
00155 // functions inherited from LSCStrategy
00156 LinearOp InvLSCStrategy::getInvBQBt(const BlockedLinearOp & A,BlockPreconditionerState & state) const
00157 {
00158    return state.getInverse("invBQBtmC");
00159 }
00160 
00161 LinearOp InvLSCStrategy::getInvBHBt(const BlockedLinearOp & A,BlockPreconditionerState & state) const
00162 {
00163    return state.getInverse("invBHBtmC");
00164 }
00165 
00166 LinearOp InvLSCStrategy::getInvF(const BlockedLinearOp & A,BlockPreconditionerState & state) const
00167 {
00168    return state.getInverse("invF");
00169 }
00170 
00171 LinearOp InvLSCStrategy::getOuterStabilization(const BlockedLinearOp & A,BlockPreconditionerState & state) const
00172 {
00173    LSCPrecondState * lscState = dynamic_cast<LSCPrecondState*>(&state);
00174    TEUCHOS_ASSERT(lscState!=0);
00175    TEUCHOS_ASSERT(lscState->isInitialized())
00176 
00177    return lscState->aiD_;
00178 }
00179 
00180 LinearOp InvLSCStrategy::getInvMass(const BlockedLinearOp & A,BlockPreconditionerState & state) const
00181 {
00182    LSCPrecondState * lscState = dynamic_cast<LSCPrecondState*>(&state);
00183    TEUCHOS_ASSERT(lscState!=0);
00184    TEUCHOS_ASSERT(lscState->isInitialized())
00185 
00186    return lscState->invMass_;
00187 }
00188 
00189 LinearOp InvLSCStrategy::getHScaling(const BlockedLinearOp & A,BlockPreconditionerState & state) const
00190 {
00191    if(hScaling_!=Teuchos::null) return hScaling_;
00192    return getInvMass(A,state);
00193 }
00194 
00196 void InvLSCStrategy::initializeState(const BlockedLinearOp & A,LSCPrecondState * state) const
00197 {
00198    Teko_DEBUG_SCOPE("InvLSCStrategy::initializeState",10);
00199 
00200    const LinearOp F  = getBlock(0,0,A);
00201    const LinearOp Bt = getBlock(0,1,A);
00202    const LinearOp B  = getBlock(1,0,A);
00203    const LinearOp C  = getBlock(1,1,A);
00204 
00205    LinearOp D = B;
00206    LinearOp G = isSymmetric_ ? Bt : adjoint(D);
00207 
00208    bool isStabilized = (not isZeroOp(C));
00209 
00210    // The logic follows like this
00211    //    if there is no mass matrix available --> build from F
00212    //    if there is a mass matrix and the inverse hasn't yet been built
00213    //       --> build from the mass matrix
00214    //    otherwise, there is already an invMass_ matrix that is appropriate
00215    //       --> use that one
00216    if(massMatrix_==Teuchos::null) {
00217       Teko_DEBUG_MSG("LSC::initializeState Build Scaling <F> type \"" 
00218                    << getDiagonalName(scaleType_) << "\"" ,1);
00219       state->invMass_ = getInvDiagonalOp(F,scaleType_);
00220    }
00221    else if(state->invMass_==Teuchos::null) {
00222       Teko_DEBUG_MSG("LSC::initializeState Build Scaling <mass> type \"" 
00223                    << getDiagonalName(scaleType_) << "\"" ,1);
00224       state->invMass_ = getInvDiagonalOp(massMatrix_,scaleType_);
00225    }
00226    // else "invMass_" should be set and there is no reason to rebuild it
00227 
00228    // compute BQBt
00229    state->BQBt_ = explicitMultiply(B,state->invMass_,Bt,state->BQBt_);
00230    Teko_DEBUG_MSG("Computed BQBt",10);
00231 
00232    // if there is no H-Scaling
00233    if(wScaling_!=Teuchos::null && hScaling_==Teuchos::null) {
00234       // from W vector build H operator scaling
00235       RCP<const Thyra::VectorBase<double> > w = wScaling_->col(0);
00236       RCP<const Thyra::VectorBase<double> > iQu 
00237             = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<double> >(state->invMass_)->getDiag();
00238       RCP<Thyra::VectorBase<double> > h = Thyra::createMember(iQu->space());
00239 
00240       Thyra::put_scalar(0.0,h.ptr());
00241       Thyra::ele_wise_prod(1.0,*w,*iQu,h.ptr());
00242       hScaling_ = Teuchos::rcp(new Thyra::DefaultDiagonalLinearOp<double>(h));
00243    } 
00244 
00245    LinearOp H = hScaling_;
00246    if(H==Teuchos::null && not isSymmetric_)
00247       H = state->invMass_;
00248 
00249    // setup the scaling operator
00250    if(H==Teuchos::null)
00251       state->BHBt_ = state->BQBt_;
00252    else {
00253       RCP<Teuchos::Time> time = Teuchos::TimeMonitor::getNewTimer("InvLSCStrategy::initializeState Build BHBt");
00254       Teuchos::TimeMonitor timer(*time);
00255 
00256       // compute BHBt
00257       state->BHBt_ = explicitMultiply(D,H,G,state->BHBt_);
00258    }
00259 
00260    // if this is a stable discretization...we are done!
00261    if(not isStabilized) {
00262       state->addInverse("BQBtmC",state->BQBt_);
00263       state->addInverse("BHBtmC",state->BHBt_);
00264       state->gamma_ = 0.0;
00265       state->alpha_ = 0.0;
00266       state->aiD_ = Teuchos::null;
00267 
00268       state->setInitialized(true);
00269 
00270       return;
00271    }
00272 
00273    // for Epetra_CrsMatrix...zero out certain rows: this ensures spectral radius is correct
00274    LinearOp modF = F;
00275    const RCP<const Epetra_Operator> epF = Thyra::get_Epetra_Operator(*F);
00276    if(epF!=Teuchos::null && rowZeroingNeeded_) {
00277       // try to get a CRS matrix
00278       const RCP<const Epetra_CrsMatrix> crsF = rcp_dynamic_cast<const Epetra_CrsMatrix>(epF);
00279 
00280       // if it is a CRS matrix get rows that need to be zeroed
00281       if(crsF!=Teuchos::null) {
00282          std::vector<int> zeroIndices;
00283           
00284          // get rows in need of zeroing
00285          Teko::Epetra::identityRowIndices(crsF->RowMap(), *crsF,zeroIndices);
00286 
00287          // build an operator that zeros those rows
00288          modF = Thyra::epetraLinearOp(rcp(new Teko::Epetra::ZeroedOperator(zeroIndices,crsF)));
00289       }
00290    }
00291 
00292    // compute gamma
00293    Teko_DEBUG_MSG("Calculating gamma",10);
00294    LinearOp iQuF = multiply(state->invMass_,modF);
00295 
00296    // do 6 power iterations to compute spectral radius: EHSST2007 Eq. 4.28
00297    Teko::LinearOp stabMatrix; // this is the pressure stabilization matrix to use
00298    state->gamma_ = std::fabs(Teko::computeSpectralRad(iQuF,5e-2,false,eigSolveParam_))/3.0; 
00299    Teko_DEBUG_MSG("Calculated gamma",10);
00300    if(userPresStabMat_!=Teuchos::null) {
00301       Teko::LinearOp invDGl = Teko::getInvDiagonalOp(userPresStabMat_);
00302       Teko::LinearOp gammaOp = multiply(invDGl,C);
00303       state->gamma_ *= std::fabs(Teko::computeSpectralRad(gammaOp,5e-2,false,eigSolveParam_));
00304       stabMatrix = userPresStabMat_;
00305    } else 
00306       stabMatrix = C;
00307 
00308    // compute alpha scaled inv(D): EHSST2007 Eq. 4.29
00309    // construct B_idF_Bt and save it for refilling later: This could reuse BQBt graph
00310    LinearOp invDiagF = getInvDiagonalOp(F);
00311    Teko::ModifiableLinearOp modB_idF_Bt = state->getInverse("BidFBt");
00312    modB_idF_Bt = explicitMultiply(B,invDiagF,Bt,modB_idF_Bt);
00313    state->addInverse("BidFBt",modB_idF_Bt);
00314    const LinearOp B_idF_Bt = modB_idF_Bt;
00315 
00316    MultiVector vec_D = getDiagonal(B_idF_Bt); // this memory could be reused
00317    update(-1.0,getDiagonal(C),1.0,vec_D); // vec_D = diag(B*inv(diag(F))*Bt)-diag(C)
00318    const LinearOp invD = buildInvDiagonal(vec_D,"inv(D)");
00319 
00320    Teko_DEBUG_MSG("Calculating alpha",10);
00321    const LinearOp BidFBtidD = multiply<double>(B_idF_Bt,invD);
00322    double num = std::fabs(Teko::computeSpectralRad(BidFBtidD,5e-2,false,eigSolveParam_));
00323    Teko_DEBUG_MSG("Calculated alpha",10);
00324    state->alpha_ = 1.0/num;
00325    state->aiD_ = Thyra::scale(state->alpha_,invD);
00326 
00327    // now build B*Q*Bt-gamma*C
00328    Teko::ModifiableLinearOp BQBtmC = state->getInverse("BQBtmC");
00329    BQBtmC = explicitAdd(state->BQBt_,scale(-state->gamma_,stabMatrix),BQBtmC);
00330    state->addInverse("BQBtmC",BQBtmC);
00331 
00332    // now build B*H*Bt-gamma*C
00333    Teko::ModifiableLinearOp BHBtmC = state->getInverse("BHBtmC");
00334    if(H==Teuchos::null)
00335       BHBtmC = BQBtmC;
00336    else {
00337       BHBtmC = explicitAdd(state->BHBt_,scale(-state->gamma_,stabMatrix),BHBtmC);
00338    }
00339    state->addInverse("BHBtmC",BHBtmC);
00340 
00341    Teko_DEBUG_MSG_BEGIN(5)
00342       DEBUG_STREAM << "LSC Gamma Parameter = " << state->gamma_ << std::endl;
00343       DEBUG_STREAM << "LSC Alpha Parameter = " << state->alpha_ << std::endl;
00344    Teko_DEBUG_MSG_END()
00345 
00346    state->setInitialized(true);
00347 }
00348 
00354 void InvLSCStrategy::computeInverses(const BlockedLinearOp & A,LSCPrecondState * state) const
00355 {
00356    Teko_DEBUG_SCOPE("InvLSCStrategy::computeInverses",10);
00357    Teko_DEBUG_EXPR(Teuchos::Time invTimer(""));
00358 
00359    const LinearOp F  = getBlock(0,0,A);
00360 
00362 
00363    // (re)build the inverse of F
00364    Teko_DEBUG_MSG("LSC::computeInverses Building inv(F)",1);
00365    Teko_DEBUG_EXPR(invTimer.start(true));
00366    InverseLinearOp invF = state->getInverse("invF");
00367    if(invF==Teuchos::null) {
00368       invF = buildInverse(*invFactoryF_,F);
00369       state->addInverse("invF",invF); 
00370    } else {
00371       rebuildInverse(*invFactoryF_,F,invF);
00372    }
00373    Teko_DEBUG_EXPR(invTimer.stop());
00374    Teko_DEBUG_MSG("LSC::computeInverses GetInvF = " << invTimer.totalElapsedTime(),1);
00375 
00377 
00378    // (re)build the inverse of BQBt 
00379    Teko_DEBUG_MSG("LSC::computeInverses Building inv(BQBtmC)",1);
00380    Teko_DEBUG_EXPR(invTimer.start(true));
00381    const LinearOp BQBt = state->getInverse("BQBtmC");
00382    InverseLinearOp invBQBt = state->getInverse("invBQBtmC");
00383    if(invBQBt==Teuchos::null) {
00384       invBQBt = buildInverse(*invFactoryS_,BQBt);
00385       state->addInverse("invBQBtmC",invBQBt); 
00386    } else {
00387       rebuildInverse(*invFactoryS_,BQBt,invBQBt);
00388    }
00389    Teko_DEBUG_EXPR(invTimer.stop());
00390    Teko_DEBUG_MSG("LSC::computeInverses GetInvBQBt = " << invTimer.totalElapsedTime(),1);
00391 
00393 
00394    // Compute the inverse of BHBt or just use BQBt
00395    ModifiableLinearOp invBHBt = state->getInverse("invBHBtmC");
00396    if(hScaling_!=Teuchos::null || not isSymmetric_) {
00397       // (re)build the inverse of BHBt 
00398       Teko_DEBUG_MSG("LSC::computeInverses Building inv(BHBtmC)",1);
00399       Teko_DEBUG_EXPR(invTimer.start(true));
00400       const LinearOp BHBt = state->getInverse("BHBtmC");
00401       if(invBHBt==Teuchos::null) {
00402          invBHBt = buildInverse(*invFactoryS_,BHBt);
00403          state->addInverse("invBHBtmC",invBHBt); 
00404       } else {
00405          rebuildInverse(*invFactoryS_,BHBt,invBHBt);
00406       }
00407       Teko_DEBUG_EXPR(invTimer.stop());
00408       Teko_DEBUG_MSG("LSC::computeInverses GetInvBHBt = " << invTimer.totalElapsedTime(),1);
00409    } 
00410    else if(invBHBt==Teuchos::null) {
00411       // just use the Q version
00412       state->addInverse("invBHBtmC",invBQBt); 
00413    }
00414 }
00415 
00417 void InvLSCStrategy::initializeFromParameterList(const Teuchos::ParameterList & pl,const InverseLibrary & invLib) 
00418 {
00419    // get string specifying inverse
00420    std::string invStr="", invVStr="", invPStr="";
00421    bool rowZeroing = true;
00422    bool useLDU = false;
00423    scaleType_ = Diagonal;
00424 
00425    // "parse" the parameter list
00426    if(pl.isParameter("Inverse Type"))
00427       invStr = pl.get<std::string>("Inverse Type");
00428    if(pl.isParameter("Inverse Velocity Type"))
00429       invVStr = pl.get<std::string>("Inverse Velocity Type");
00430    if(pl.isParameter("Inverse Pressure Type")) 
00431       invPStr = pl.get<std::string>("Inverse Pressure Type");
00432    if(pl.isParameter("Ignore Boundary Rows"))
00433       rowZeroing = pl.get<bool>("Ignore Boundary Rows");
00434    if(pl.isParameter("Use LDU"))
00435       useLDU = pl.get<bool>("Use LDU");
00436    if(pl.isParameter("Use Mass Scaling"))
00437       useMass_ = pl.get<bool>("Use Mass Scaling");
00438    // if(pl.isParameter("Use Lumping"))
00439    //    useLumping_ = pl.get<bool>("Use Lumping");
00440    if(pl.isParameter("Use W-Scaling"))
00441       useWScaling_ = pl.get<bool>("Use W-Scaling");
00442    if(pl.isParameter("Eigen Solver Iterations"))
00443       eigSolveParam_ = pl.get<int>("Eigen Solver Iterations");
00444    if(pl.isParameter("Scaling Type")) {
00445       scaleType_ = getDiagonalType(pl.get<std::string>("Scaling Type"));
00446       TEUCHOS_TEST_FOR_EXCEPT(scaleType_==NotDiag);
00447    }
00448 
00449    Teko_DEBUG_MSG_BEGIN(5)
00450       DEBUG_STREAM << "LSC Inverse Strategy Parameters: " << std::endl;
00451       DEBUG_STREAM << "   inv type   = \"" << invStr  << "\"" << std::endl;
00452       DEBUG_STREAM << "   inv v type = \"" << invVStr << "\"" << std::endl;
00453       DEBUG_STREAM << "   inv p type = \"" << invPStr << "\"" << std::endl;
00454       DEBUG_STREAM << "   bndry rows = " << rowZeroing << std::endl;
00455       DEBUG_STREAM << "   use ldu    = " << useLDU << std::endl;
00456       DEBUG_STREAM << "   use mass    = " << useMass_ << std::endl;
00457       DEBUG_STREAM << "   use w-scaling    = " << useWScaling_ << std::endl;
00458       DEBUG_STREAM << "   scale type    = " << getDiagonalName(scaleType_) << std::endl;
00459       DEBUG_STREAM << "LSC  Inverse Strategy Parameter list: " << std::endl;
00460       pl.print(DEBUG_STREAM);
00461    Teko_DEBUG_MSG_END()
00462 
00463    // set defaults as needed
00464    if(invStr=="") invStr = "Amesos";
00465    if(invVStr=="") invVStr = invStr;
00466    if(invPStr=="") invPStr = invStr;
00467 
00468    // build velocity inverse factory
00469    invFactoryF_ = invLib.getInverseFactory(invVStr);
00470    invFactoryS_ = invFactoryF_; // by default these are the same
00471    if(invVStr!=invPStr) // if different, build pressure inverse factory
00472       invFactoryS_ = invLib.getInverseFactory(invPStr);
00473 
00474    // set other parameters
00475    setUseFullLDU(useLDU);
00476    setRowZeroing(rowZeroing);
00477 
00478    if(useMass_) {
00479       Teuchos::RCP<Teko::RequestHandler> rh = getRequestHandler();
00480       rh->preRequest<Teko::LinearOp>(Teko::RequestMesg("Velocity Mass Matrix"));
00481       Teko::LinearOp mass 
00482             = rh->request<Teko::LinearOp>(Teko::RequestMesg("Velocity Mass Matrix"));
00483       setMassMatrix(mass);
00484    }
00485 
00486 }
00487 
00489 Teuchos::RCP<Teuchos::ParameterList> InvLSCStrategy::getRequestedParameters() const 
00490 {
00491    Teuchos::RCP<Teuchos::ParameterList> result;
00492    Teuchos::RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList());
00493 
00494    // grab parameters from F solver
00495    RCP<Teuchos::ParameterList> fList = invFactoryF_->getRequestedParameters();
00496    if(fList!=Teuchos::null) {
00497       Teuchos::ParameterList::ConstIterator itr;
00498       for(itr=fList->begin();itr!=fList->end();++itr)
00499          pl->setEntry(itr->first,itr->second);
00500       result = pl;
00501    }
00502 
00503    // grab parameters from S solver
00504    RCP<Teuchos::ParameterList> sList = invFactoryS_->getRequestedParameters();
00505    if(sList!=Teuchos::null) {
00506       Teuchos::ParameterList::ConstIterator itr;
00507       for(itr=sList->begin();itr!=sList->end();++itr)
00508          pl->setEntry(itr->first,itr->second);
00509       result = pl;
00510    }
00511 
00512    // use the mass matrix
00513    if(useWScaling_) {
00514       pl->set<Teko::LinearOp>("W-Scaling Vector", Teuchos::null,"W-Scaling Vector");
00515       result = pl;
00516    }
00517 
00518    return result;
00519 }
00520 
00522 bool InvLSCStrategy::updateRequestedParameters(const Teuchos::ParameterList & pl) 
00523 {
00524    Teko_DEBUG_SCOPE("InvLSCStrategy::updateRequestedParameters",10);
00525    bool result = true;
00526  
00527    // update requested parameters in solvers
00528    result &= invFactoryF_->updateRequestedParameters(pl);
00529    result &= invFactoryS_->updateRequestedParameters(pl);
00530 
00531    // use W scaling matrix
00532    if(useWScaling_) {
00533       Teko::MultiVector wScale = pl.get<Teko::MultiVector>("W-Scaling Vector");
00534 
00535       if(wScale==Teuchos::null)
00536          result &= false;
00537       else
00538          setWScaling(wScale);
00539    }
00540 
00541    return result;
00542 }
00543 
00544 } // end namespace NS
00545 } // end namespace Teko
 All Classes Files Functions Variables