Teko Version of the Day
Teko_SIMPLEPreconditionerFactory.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 "Teko_SIMPLEPreconditionerFactory.hpp"
00048 
00049 #include "Teko_Utilities.hpp"
00050 #include "Teko_InverseFactory.hpp"
00051 #include "Teko_BlockLowerTriInverseOp.hpp"
00052 #include "Teko_BlockUpperTriInverseOp.hpp"
00053 #include "Teko_DiagonalPreconditionerFactory.hpp"
00054 
00055 #include "Teuchos_Time.hpp"
00056 #include "Epetra_FECrsGraph.h"
00057 #include "Epetra_FECrsMatrix.h"
00058 #include "EpetraExt_PointToBlockDiagPermute.h"
00059 #include "EpetraExt_MatrixMatrix.h"
00060 #include "Thyra_EpetraOperatorWrapper.hpp"
00061 #include "Thyra_EpetraLinearOp.hpp"
00062 
00063 
00064 using Teuchos::RCP;
00065 
00066 namespace Teko {
00067 namespace NS {
00068 
00069 // Constructor definition
00070 SIMPLEPreconditionerFactory
00071    ::SIMPLEPreconditionerFactory(const RCP<InverseFactory> & inverse,
00072                                  double alpha)
00073    : invVelFactory_(inverse), invPrsFactory_(inverse), alpha_(alpha), fInverseType_(Diagonal), useMass_(false)
00074 { }
00075 
00076 SIMPLEPreconditionerFactory
00077    ::SIMPLEPreconditionerFactory(const RCP<InverseFactory> & invVFact,
00078                                  const RCP<InverseFactory> & invPFact,
00079                                  double alpha)
00080    : invVelFactory_(invVFact), invPrsFactory_(invPFact), alpha_(alpha), fInverseType_(Diagonal), useMass_(false)
00081 { }
00082 
00083 SIMPLEPreconditionerFactory::SIMPLEPreconditionerFactory()
00084    : alpha_(1.0), fInverseType_(Diagonal), useMass_(false)
00085 { }
00086 
00087 // Use the factory to build the preconditioner (this is where the work goes)
00088 LinearOp SIMPLEPreconditionerFactory
00089    ::buildPreconditionerOperator(BlockedLinearOp & blockOp,
00090                                  BlockPreconditionerState & state) const
00091 {
00092    Teko_DEBUG_SCOPE("SIMPLEPreconditionerFactory::buildPreconditionerOperator",10);
00093    Teko_DEBUG_EXPR(Teuchos::Time timer(""));
00094 
00095    int rows = blockRowCount(blockOp);
00096    int cols = blockColCount(blockOp);
00097  
00098    TEUCHOS_ASSERT(rows==2); // sanity checks
00099    TEUCHOS_ASSERT(cols==2);
00100 
00101    bool buildExplicitSchurComplement = true;
00102 
00103    // extract subblocks
00104    const LinearOp F  = getBlock(0,0,blockOp);
00105    const LinearOp Bt = getBlock(0,1,blockOp);
00106    const LinearOp B  = getBlock(1,0,blockOp);
00107    const LinearOp C  = getBlock(1,1,blockOp);
00108 
00109    LinearOp matF = F;
00110    if(useMass_) {
00111       TEUCHOS_ASSERT(massMatrix_!=Teuchos::null);
00112       matF = massMatrix_;
00113    }
00114 
00115    // get approximation of inv(F) name H
00116    std::string fApproxStr = "<error>";
00117    LinearOp H;
00118    if(fInverseType_==NotDiag) {
00119       H = buildInverse(*customHFactory_,matF);
00120       fApproxStr = customHFactory_->toString();
00121 
00122       // since H is now implicit, we must build an implicit Schur complement
00123       buildExplicitSchurComplement = false;
00124    }
00125    else if(fInverseType_==BlkDiag) {
00126      // Block diagonal approximation for H
00127      DiagonalPreconditionerFactory Hfact;
00128      DiagonalPrecondState Hstate;
00129      Hfact.initializeFromParameterList(BlkDiagList_);           
00130      H = Hfact.buildPreconditionerOperator(matF,Hstate); 
00131 
00132 /*
00133      // Get a FECrsMarix out of the BDP
00134      RCP<Epetra_FECrsMatrix> Hcrs=rcp(Hstate.BDP_->CreateFECrsMatrix());
00135      H=Thyra::epetraLinearOp(Hcrs);
00136 */
00137 
00138      buildExplicitSchurComplement = true; // NTS: Do I need this? 
00139                                           // Answer - no, but it is documenting whats going on here.
00140    }
00141    else {
00142       // get generic diagonal
00143       H = getInvDiagonalOp(matF,fInverseType_);
00144       fApproxStr = getDiagonalName(fInverseType_);
00145    }
00146 
00147    // adjust H for time scaling if it is a mass matrix
00148    if(useMass_) {
00149       RCP<const Teuchos::ParameterList> pl = state.getParameterList();
00150 
00151       if(pl->isParameter("stepsize")) {
00152          // get the step size
00153          double stepsize = pl->get<double>("stepsize");
00154 
00155          // scale by stepsize only if it is larger than 0
00156          if(stepsize>0.0)
00157             H = scale(stepsize,H);
00158       }
00159    }
00160 
00161    // build approximate Schur complement: hatS = -C + B*H*Bt
00162    LinearOp HBt, hatS;
00163 
00164    if(buildExplicitSchurComplement) {
00165       ModifiableLinearOp & mHBt = state.getModifiableOp("HBt");
00166       ModifiableLinearOp & mhatS = state.getModifiableOp("hatS");
00167       ModifiableLinearOp & BHBt = state.getModifiableOp("BHBt");
00168 
00169       // build H*Bt
00170       mHBt = explicitMultiply(H,Bt,mHBt);
00171       HBt = mHBt;
00172 
00173       // build B*H*Bt
00174       BHBt = explicitMultiply(B,HBt,BHBt);
00175 
00176       // build C-B*H*Bt
00177       mhatS = explicitAdd(C,scale(-1.0,BHBt),mhatS);
00178       hatS = mhatS;
00179    }
00180    else {
00181       // build an implicit Schur complement
00182       HBt = multiply(H,Bt);
00183       
00184       hatS = add(C,scale(-1.0,multiply(B,HBt)));
00185    }
00186 
00187    // build the inverse for F 
00188    ModifiableLinearOp & invF = state.getModifiableOp("invF");
00189    if(invF==Teuchos::null)
00190       invF = buildInverse(*invVelFactory_,F);
00191    else 
00192       rebuildInverse(*invVelFactory_,F,invF);
00193 
00194    // build the approximate Schur complement: This is inefficient! FIXME
00195    ModifiableLinearOp & invS = state.getModifiableOp("invS");
00196    if(invS==Teuchos::null)
00197       invS = buildInverse(*invPrsFactory_,hatS);
00198    else
00199       rebuildInverse(*invPrsFactory_,hatS,invS);
00200 
00201    std::vector<LinearOp> invDiag(2); // vector storing inverses
00202 
00203    // build lower triangular inverse matrix
00204    BlockedLinearOp L = zeroBlockedOp(blockOp);
00205    setBlock(1,0,L,B);
00206    endBlockFill(L);
00207 
00208    invDiag[0] = invF;
00209    invDiag[1] = invS;
00210    LinearOp invL = createBlockLowerTriInverseOp(L,invDiag);
00211 
00212    // build upper triangular matrix
00213    BlockedLinearOp U = zeroBlockedOp(blockOp);
00214    setBlock(0,1,U,scale(1.0/alpha_,HBt));
00215    endBlockFill(U);
00216 
00217    invDiag[0] = identity(rangeSpace(invF));
00218    invDiag[1] = scale(alpha_,identity(rangeSpace(invS)));
00219    LinearOp invU = createBlockUpperTriInverseOp(U,invDiag);
00220 
00221    // return implicit product operator
00222    return multiply(invU,invL,"SIMPLE_"+fApproxStr);
00223 }
00224 
00226 void SIMPLEPreconditionerFactory::initializeFromParameterList(const Teuchos::ParameterList & pl)
00227 {
00228    RCP<const InverseLibrary> invLib = getInverseLibrary();
00229 
00230    // default conditions
00231    useMass_ = false;
00232    customHFactory_ = Teuchos::null;
00233    fInverseType_ = Diagonal;
00234   
00235    // get string specifying inverse
00236    std::string invStr="", invVStr="", invPStr="";
00237    alpha_ = 1.0;
00238 
00239    // "parse" the parameter list
00240    if(pl.isParameter("Inverse Type"))
00241       invStr = pl.get<std::string>("Inverse Type");
00242    if(pl.isParameter("Inverse Velocity Type"))
00243      invVStr = pl.get<std::string>("Inverse Velocity Type");
00244    if(pl.isParameter("Inverse Pressure Type")) 
00245      invPStr = pl.get<std::string>("Inverse Pressure Type");
00246    if(pl.isParameter("Alpha"))
00247      alpha_ = pl.get<double>("Alpha");
00248    if(pl.isParameter("Explicit Velocity Inverse Type")) {
00249       std::string fInverseStr = pl.get<std::string>("Explicit Velocity Inverse Type");
00250 
00251       // build inverse types
00252       fInverseType_ = getDiagonalType(fInverseStr);
00253       if(fInverseType_==NotDiag)
00254          customHFactory_ = invLib->getInverseFactory(fInverseStr);
00255 
00256       // Grab the sublist if we're using the block diagonal
00257       if(fInverseType_==BlkDiag)
00258   BlkDiagList_=pl.sublist("H options");      
00259    }
00260    if(pl.isParameter("Use Mass Scaling"))
00261       useMass_ = pl.get<bool>("Use Mass Scaling");
00262 
00263    Teko_DEBUG_MSG_BEGIN(5)
00264       DEBUG_STREAM << "SIMPLE Parameters: " << std::endl;
00265       DEBUG_STREAM << "   inv type    = \"" << invStr  << "\"" << std::endl;
00266       DEBUG_STREAM << "   inv v type  = \"" << invVStr << "\"" << std::endl;
00267       DEBUG_STREAM << "   inv p type  = \"" << invPStr << "\"" << std::endl;
00268       DEBUG_STREAM << "   alpha       = " << alpha_ << std::endl;
00269       DEBUG_STREAM << "   use mass    = " << useMass_ << std::endl;
00270       DEBUG_STREAM << "   vel scaling = " << getDiagonalName(fInverseType_) << std::endl;
00271       DEBUG_STREAM << "SIMPLE Parameter list: " << std::endl;
00272       pl.print(DEBUG_STREAM);
00273    Teko_DEBUG_MSG_END()
00274 
00275    // set defaults as needed
00276    if(invStr=="") invStr = "Amesos";
00277    if(invVStr=="") invVStr = invStr;
00278    if(invPStr=="") invPStr = invStr;
00279 
00280    //  two inverse factory objects
00281    RCP<InverseFactory> invVFact, invPFact;
00282 
00283    // build velocity inverse factory
00284    invVFact = invLib->getInverseFactory(invVStr);
00285    invPFact = invVFact; // by default these are the same
00286    if(invVStr!=invPStr) // if different, build pressure inverse factory
00287       invPFact = invLib->getInverseFactory(invPStr);
00288 
00289    // based on parameter type build a strategy
00290    invVelFactory_ = invVFact; 
00291    invPrsFactory_ = invPFact;
00292 
00293    if(useMass_) {
00294       Teuchos::RCP<Teko::RequestHandler> rh = getRequestHandler();
00295       rh->preRequest<Teko::LinearOp>(Teko::RequestMesg("Velocity Mass Matrix"));
00296       Teko::LinearOp mass 
00297             = rh->request<Teko::LinearOp>(Teko::RequestMesg("Velocity Mass Matrix"));
00298       setMassMatrix(mass);
00299    }
00300 }
00301 
00303 Teuchos::RCP<Teuchos::ParameterList> SIMPLEPreconditionerFactory::getRequestedParameters() const 
00304 {
00305    Teuchos::RCP<Teuchos::ParameterList> result;
00306    Teuchos::RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList());
00307 
00308    // grab parameters from F solver
00309    RCP<Teuchos::ParameterList> vList = invVelFactory_->getRequestedParameters();
00310    if(vList!=Teuchos::null) {
00311       Teuchos::ParameterList::ConstIterator itr;
00312       for(itr=vList->begin();itr!=vList->end();++itr)
00313          pl->setEntry(itr->first,itr->second);
00314       result = pl;
00315    }
00316 
00317    // grab parameters from S solver
00318    RCP<Teuchos::ParameterList> pList = invPrsFactory_->getRequestedParameters();
00319    if(pList!=Teuchos::null) {
00320       Teuchos::ParameterList::ConstIterator itr;
00321       for(itr=pList->begin();itr!=pList->end();++itr)
00322          pl->setEntry(itr->first,itr->second);
00323       result = pl;
00324    }
00325 
00326    // grab parameters from S solver
00327    if(customHFactory_!=Teuchos::null) {
00328       RCP<Teuchos::ParameterList> hList = customHFactory_->getRequestedParameters();
00329       if(hList!=Teuchos::null) {
00330          Teuchos::ParameterList::ConstIterator itr;
00331          for(itr=hList->begin();itr!=hList->end();++itr)
00332             pl->setEntry(itr->first,itr->second);
00333          result = pl;
00334       }
00335    }
00336 
00337    return result;
00338 }
00339 
00341 bool SIMPLEPreconditionerFactory::updateRequestedParameters(const Teuchos::ParameterList & pl) 
00342 {
00343    Teko_DEBUG_SCOPE("InvLSCStrategy::updateRequestedParameters",10);
00344    bool result = true;
00345  
00346    // update requested parameters in solvers
00347    result &= invVelFactory_->updateRequestedParameters(pl);
00348    result &= invPrsFactory_->updateRequestedParameters(pl);
00349    if(customHFactory_!=Teuchos::null)
00350       result &= customHFactory_->updateRequestedParameters(pl);
00351 
00352    return result;
00353 }
00354 
00355 } // end namespace NS
00356 } // end namespace Teko
 All Classes Files Functions Variables