Teko Version of the Day
Teko_PresLaplaceLSCStrategy.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_PresLaplaceLSCStrategy.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 
00054 #include "Teuchos_Time.hpp"
00055 #include "Teuchos_TimeMonitor.hpp"
00056 
00057 // Teko includes
00058 #include "Teko_Utilities.hpp"
00059 #include "NS/Teko_LSCPreconditionerFactory.hpp"
00060 
00061 using Teuchos::RCP;
00062 using Teuchos::rcp_dynamic_cast;
00063 using Teuchos::rcp_const_cast;
00064 
00065 namespace Teko {
00066 namespace NS {
00067 
00069 // PresLaplaceLSCStrategy Implementation
00071 
00072 // constructors
00074 PresLaplaceLSCStrategy::PresLaplaceLSCStrategy()
00075    : invFactoryV_(Teuchos::null), invFactoryP_(Teuchos::null)
00076    , eigSolveParam_(5), useFullLDU_(false), useMass_(false), scaleType_(AbsRowSum)
00077 { }
00078 
00079 PresLaplaceLSCStrategy::PresLaplaceLSCStrategy(const Teuchos::RCP<InverseFactory> & factory)
00080    : invFactoryV_(factory), invFactoryP_(factory)
00081    , eigSolveParam_(5), useFullLDU_(false), useMass_(false), scaleType_(AbsRowSum)
00082 { }
00083 
00084 PresLaplaceLSCStrategy::PresLaplaceLSCStrategy(const Teuchos::RCP<InverseFactory> & invFactF,
00085                                                const Teuchos::RCP<InverseFactory> & invFactS)
00086    : invFactoryV_(invFactF), invFactoryP_(invFactS)
00087    , eigSolveParam_(5), useFullLDU_(false), useMass_(false), scaleType_(AbsRowSum)
00088 { }
00089 
00091 
00092 void PresLaplaceLSCStrategy::buildState(BlockedLinearOp & A,BlockPreconditionerState & state) const
00093 {
00094    Teko_DEBUG_SCOPE("PresLaplaceLSCStrategy::buildState",10);
00095 
00096    LSCPrecondState * lscState = dynamic_cast<LSCPrecondState*>(&state);
00097    TEUCHOS_ASSERT(lscState!=0);
00098 
00099    // if neccessary save state information
00100    if(not lscState->isInitialized()) {
00101       Teko_DEBUG_EXPR(Teuchos::Time timer(""));
00102 
00103       // construct operators
00104       {
00105          Teko_DEBUG_SCOPE("PL-LSC::buildState constructing operators",1);
00106          Teko_DEBUG_EXPR(timer.start(true));
00107 
00108          initializeState(A,lscState);
00109 
00110          Teko_DEBUG_EXPR(timer.stop());
00111          Teko_DEBUG_MSG("PL-LSC::buildState BuildOpsTime = " << timer.totalElapsedTime(),1);
00112       }
00113 
00114       // Build the inverses
00115       {
00116          Teko_DEBUG_SCOPE("PL-LSC::buildState calculating inverses",1);
00117          Teko_DEBUG_EXPR(timer.start(true));
00118 
00119          computeInverses(A,lscState);
00120 
00121          Teko_DEBUG_EXPR(timer.stop());
00122          Teko_DEBUG_MSG("PL-LSC::buildState BuildInvTime = " << timer.totalElapsedTime(),1);
00123       }
00124    }
00125 }
00126 
00127 // functions inherited from LSCStrategy
00128 LinearOp PresLaplaceLSCStrategy::getInvBQBt(const BlockedLinearOp & A,BlockPreconditionerState & state) const
00129 {
00130    return state.getModifiableOp("invPresLap");
00131 }
00132 
00133 LinearOp PresLaplaceLSCStrategy::getInvBHBt(const BlockedLinearOp & A,BlockPreconditionerState & state) const
00134 {
00135    return state.getModifiableOp("invPresLap");
00136 }
00137 
00138 LinearOp PresLaplaceLSCStrategy::getInvF(const BlockedLinearOp & A,BlockPreconditionerState & state) const
00139 {
00140    return state.getModifiableOp("invF");
00141 }
00142 
00143 // LinearOp PresLaplaceLSCStrategy::getInvAlphaD(const BlockedLinearOp & A,BlockPreconditionerState & state) const
00144 LinearOp PresLaplaceLSCStrategy::getOuterStabilization(const BlockedLinearOp & A,BlockPreconditionerState & state) const
00145 {
00146    LSCPrecondState * lscState = dynamic_cast<LSCPrecondState*>(&state);
00147    TEUCHOS_ASSERT(lscState!=0);
00148    TEUCHOS_ASSERT(lscState->isInitialized())
00149 
00150    return lscState->aiD_;
00151 }
00152 
00153 LinearOp PresLaplaceLSCStrategy::getInvMass(const BlockedLinearOp & A,BlockPreconditionerState & state) const
00154 {
00155    LSCPrecondState * lscState = dynamic_cast<LSCPrecondState*>(&state);
00156    TEUCHOS_ASSERT(lscState!=0);
00157    TEUCHOS_ASSERT(lscState->isInitialized())
00158 
00159    return lscState->invMass_;
00160 }
00161 
00162 LinearOp PresLaplaceLSCStrategy::getHScaling(const BlockedLinearOp & A,BlockPreconditionerState & state) const
00163 {
00164    return getInvMass(A,state);
00165 }
00166 
00168 void PresLaplaceLSCStrategy::initializeState(const BlockedLinearOp & A,LSCPrecondState * state) const
00169 {
00170    Teko_DEBUG_SCOPE("PresLaplaceLSCStrategy::initializeState",10);
00171 
00172    std::string velMassStr = getVelocityMassString();
00173 
00174    const LinearOp F  = getBlock(0,0,A);
00175    const LinearOp Bt = getBlock(0,1,A);
00176    const LinearOp B  = getBlock(1,0,A);
00177    const LinearOp C  = getBlock(1,1,A);
00178 
00179    LinearOp D = B;
00180    LinearOp G = Bt;
00181 
00182    bool isStabilized = (not isZeroOp(C));
00183 
00184    // grab operators from state object
00185    LinearOp massMatrix = state->getLinearOp(velMassStr);
00186 
00187    // The logic follows like this
00188    //    if there is no mass matrix available --> build from F
00189    //    if there is a mass matrix and the inverse hasn't yet been built
00190    //       --> build from the mass matrix
00191    //    otherwise, there is already an invMass_ matrix that is appropriate
00192    //       --> use that one
00193    if(massMatrix==Teuchos::null) {
00194       Teko_DEBUG_MSG("PL-LSC::initializeState Build Scaling <F> type \"" 
00195                    << getDiagonalName(scaleType_) << "\"" ,1);
00196       state->invMass_ = getInvDiagonalOp(F,scaleType_);
00197    }
00198    else if(state->invMass_==Teuchos::null) {
00199       Teko_DEBUG_MSG("PL-LSC::initializeState Build Scaling <mass> type \"" 
00200                    << getDiagonalName(scaleType_) << "\"" ,1);
00201       state->invMass_ = getInvDiagonalOp(massMatrix,scaleType_);
00202    }
00203    // else "invMass_" should be set and there is no reason to rebuild it
00204 
00205    // if this is a stable discretization...we are done!
00206    if(not isStabilized) {
00207       state->aiD_ = Teuchos::null;
00208 
00209       state->setInitialized(true);
00210 
00211       return;
00212    }
00213 
00214    // compute alpha scaled inv(D): EHSST2007 Eq. 4.29
00215    // construct B_idF_Bt and save it for refilling later: This could reuse BQBt graph
00216    LinearOp invDiagF = getInvDiagonalOp(F);
00217    Teko::ModifiableLinearOp & modB_idF_Bt = state->getModifiableOp("BidFBt");
00218    modB_idF_Bt = explicitMultiply(B,invDiagF,Bt,modB_idF_Bt);
00219    const LinearOp B_idF_Bt = modB_idF_Bt;
00220 
00221    MultiVector vec_D = getDiagonal(B_idF_Bt); // this memory could be reused
00222    update(-1.0,getDiagonal(C),1.0,vec_D); // vec_D = diag(B*inv(diag(F))*Bt)-diag(C)
00223    const LinearOp invD = buildInvDiagonal(vec_D,"inv(D)");
00224 
00225    Teko_DEBUG_MSG("Calculating alpha",10);
00226    const LinearOp BidFBtidD = multiply<double>(B_idF_Bt,invD);
00227    double num = std::fabs(Teko::computeSpectralRad(BidFBtidD,5e-2,false,eigSolveParam_));
00228    Teko_DEBUG_MSG("Calculated alpha",10);
00229    state->alpha_ = 1.0/num;
00230    state->aiD_ = Thyra::scale(state->alpha_,invD);
00231 
00232    Teko_DEBUG_MSG_BEGIN(5)
00233       DEBUG_STREAM << "PL-LSC Alpha Parameter = " << state->alpha_ << std::endl;
00234    Teko_DEBUG_MSG_END()
00235 
00236    state->setInitialized(true);
00237 }
00238 
00244 void PresLaplaceLSCStrategy::computeInverses(const BlockedLinearOp & A,LSCPrecondState * state) const
00245 {
00246    Teko_DEBUG_SCOPE("PresLaplaceLSCStrategy::computeInverses",10);
00247    Teko_DEBUG_EXPR(Teuchos::Time invTimer(""));
00248 
00249    std::string presLapStr  = getPressureLaplaceString();
00250 
00251    const LinearOp F  = getBlock(0,0,A);
00252    const LinearOp presLap = state->getLinearOp(presLapStr);
00253 
00255 
00256    // (re)build the inverse of F
00257    Teko_DEBUG_MSG("PL-LSC::computeInverses Building inv(F)",1);
00258    Teko_DEBUG_EXPR(invTimer.start(true));
00259    ModifiableLinearOp & invF = state->getModifiableOp("invF");
00260    if(invF==Teuchos::null) {
00261       invF = buildInverse(*invFactoryV_,F);
00262    } else {
00263       rebuildInverse(*invFactoryV_,F,invF);
00264    }
00265    Teko_DEBUG_EXPR(invTimer.stop());
00266    Teko_DEBUG_MSG("PL-LSC::computeInverses GetInvF = " << invTimer.totalElapsedTime(),1);
00267 
00269 
00270    // (re)build the inverse of P
00271    Teko_DEBUG_MSG("PL-LSC::computeInverses Building inv(PresLap)",1);
00272    Teko_DEBUG_EXPR(invTimer.start(true));
00273    ModifiableLinearOp & invPresLap = state->getModifiableOp("invPresLap");
00274    if(invPresLap==Teuchos::null) {
00275       invPresLap = buildInverse(*invFactoryP_,presLap);
00276    } else {
00277       // not need because the pressure laplacian never changes
00278       // rebuildInverse(*invFactoryP_,presLap,invPresLap);
00279    }
00280    Teko_DEBUG_EXPR(invTimer.stop());
00281    Teko_DEBUG_MSG("PL-LSC::computeInverses GetInvBQBt = " << invTimer.totalElapsedTime(),1);
00282 }
00283 
00285 void PresLaplaceLSCStrategy::initializeFromParameterList(const Teuchos::ParameterList & pl,const InverseLibrary & invLib) 
00286 {
00287    // get string specifying inverse
00288    std::string invStr="Amesos", invVStr="", invPStr="";
00289    bool useLDU = false;
00290    scaleType_ = AbsRowSum;
00291 
00292    // "parse" the parameter list
00293    if(pl.isParameter("Inverse Type"))
00294       invStr = pl.get<std::string>("Inverse Type");
00295    if(pl.isParameter("Inverse Velocity Type"))
00296       invVStr = pl.get<std::string>("Inverse Velocity Type");
00297    if(pl.isParameter("Inverse Pressure Type")) 
00298       invPStr = pl.get<std::string>("Inverse Pressure Type");
00299    if(pl.isParameter("Use LDU"))
00300       useLDU = pl.get<bool>("Use LDU");
00301    if(pl.isParameter("Use Mass Scaling"))
00302       useMass_ = pl.get<bool>("Use Mass Scaling");
00303    if(pl.isParameter("Eigen Solver Iterations"))
00304       eigSolveParam_ = pl.get<int>("Eigen Solver Iterations");
00305    if(pl.isParameter("Scaling Type")) {
00306       scaleType_ = getDiagonalType(pl.get<std::string>("Scaling Type"));
00307       TEUCHOS_TEST_FOR_EXCEPT(scaleType_==NotDiag);
00308    }
00309 
00310    // set defaults as needed
00311    if(invVStr=="") invVStr = invStr;
00312    if(invPStr=="") invPStr = invStr;
00313 
00314    Teko_DEBUG_MSG_BEGIN(5)
00315       DEBUG_STREAM << "LSC Inverse Strategy Parameters: " << std::endl;
00316       DEBUG_STREAM << "   inv v type = \"" << invVStr << "\"" << std::endl;
00317       DEBUG_STREAM << "   inv p type = \"" << invPStr << "\"" << std::endl;
00318       DEBUG_STREAM << "   use ldu    = " << useLDU << std::endl;
00319       DEBUG_STREAM << "   use mass    = " << useMass_ << std::endl;
00320       DEBUG_STREAM << "   scale type    = " << getDiagonalName(scaleType_) << std::endl;
00321       DEBUG_STREAM << "LSC  Pressure Laplace Strategy Parameter list: " << std::endl;
00322       pl.print(DEBUG_STREAM);
00323    Teko_DEBUG_MSG_END()
00324 
00325    // build velocity inverse factory
00326    invFactoryV_ = invLib.getInverseFactory(invVStr);
00327    invFactoryP_ = invFactoryV_; // by default these are the same
00328    if(invVStr!=invPStr) // if different, build pressure inverse factory
00329       invFactoryP_ = invLib.getInverseFactory(invPStr);
00330 
00331    // set other parameters
00332    setUseFullLDU(useLDU);
00333 }
00334 
00336 Teuchos::RCP<Teuchos::ParameterList> PresLaplaceLSCStrategy::getRequestedParameters() const 
00337 {
00338    Teko_DEBUG_SCOPE("PresLaplaceLSCStrategy::getRequestedParameters",10);
00339    Teuchos::RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList());
00340 
00341    // grab parameters from F solver
00342    RCP<Teuchos::ParameterList> fList = invFactoryV_->getRequestedParameters();
00343    if(fList!=Teuchos::null) {
00344       Teuchos::ParameterList::ConstIterator itr;
00345       for(itr=fList->begin();itr!=fList->end();++itr)
00346          pl->setEntry(itr->first,itr->second);
00347    }
00348 
00349    // grab parameters from S solver
00350    RCP<Teuchos::ParameterList> sList = invFactoryP_->getRequestedParameters();
00351    if(sList!=Teuchos::null) {
00352       Teuchos::ParameterList::ConstIterator itr;
00353       for(itr=sList->begin();itr!=sList->end();++itr)
00354          pl->setEntry(itr->first,itr->second);
00355    }
00356 
00357    // use the mass matrix
00358    if(useMass_)
00359       pl->set(getVelocityMassString(), false,"Velocity mass matrix");
00360    pl->set(getPressureLaplaceString(), false,"Pressure Laplacian matrix");
00361 
00362    return pl;
00363 }
00364 
00366 bool PresLaplaceLSCStrategy::updateRequestedParameters(const Teuchos::ParameterList & pl) 
00367 {
00368    Teko_DEBUG_SCOPE("PresLaplaceLSCStrategy::updateRequestedParameters",10);
00369    bool result = true;
00370  
00371    // update requested parameters in solvers
00372    result &= invFactoryV_->updateRequestedParameters(pl);
00373    result &= invFactoryP_->updateRequestedParameters(pl);
00374 
00375    Teuchos::ParameterList hackList(pl);
00376 
00377    // get required operator acknowledgment...user must set these to true
00378    bool plo = hackList.get<bool>(getPressureLaplaceString(),false);
00379  
00380    bool vmo = true;
00381    if(useMass_)
00382       vmo = hackList.get<bool>(getVelocityMassString(),false);
00383 
00384    if(not plo)  { Teko_DEBUG_MSG("User must acknowledge the use of the \"" << getPressureLaplaceString() << "\"!",0); }
00385    if(not vmo)  { Teko_DEBUG_MSG("User must acknowledge the use of the \"" << getVelocityMassString() << "\"!",0); }
00386 
00387    result &= (plo & vmo);
00388 
00389    return result;
00390 }
00391 
00392 } // end namespace NS
00393 } // end namespace Teko
 All Classes Files Functions Variables