Teko Version of the Day
Teko_PCDStrategy.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_PCDStrategy.hpp"
00048 
00049 #include "Teuchos_TimeMonitor.hpp"
00050 #include "Teko_Utilities.hpp"
00051 
00052 namespace Teko {
00053 namespace NS {
00054 
00055 using Teuchos::TimeMonitor;
00056 
00057 Teuchos::RCP<Teuchos::Time> PCDStrategy::initTimer_;
00058 Teuchos::RCP<Teuchos::Time> PCDStrategy::invSTimer_;
00059 Teuchos::RCP<Teuchos::Time> PCDStrategy::invFTimer_;
00060 Teuchos::RCP<Teuchos::Time> PCDStrategy::opsTimer_;
00061 
00062 void PCDStrategy::buildTimers()
00063 {
00064    if(initTimer_==Teuchos::null)
00065       initTimer_ = TimeMonitor::getNewTimer("PCDStrategy::initializePrec");
00066 
00067    if(invSTimer_==Teuchos::null)
00068       invSTimer_ = TimeMonitor::getNewTimer("PCDStrategy::initializePrec invS");
00069 
00070    if(invFTimer_==Teuchos::null)
00071       invFTimer_ = TimeMonitor::getNewTimer("PCDStrategy::initializePrec invF");
00072 
00073    if(opsTimer_==Teuchos::null)
00074       opsTimer_ = TimeMonitor::getNewTimer("PCDStrategy::initializePrec buildOps");
00075 }
00076 
00077 PCDStrategy::PCDStrategy() : massInverseType_(Diagonal), schurCompOrdering_(false)
00078 { 
00079    pcdParams_ = Teuchos::rcp(new Teuchos::ParameterList);
00080    lapParams_ = Teuchos::rcp(new Teuchos::ParameterList);
00081 
00082    lapParams_->set("Name",getPressureLaplaceString());
00083    pcdParams_->set("Name",getPCDString());
00084 
00085    buildTimers();
00086 }
00087 
00089 PCDStrategy::PCDStrategy(const Teuchos::RCP<InverseFactory> & invFA,
00090                                              const Teuchos::RCP<InverseFactory> & invS)
00091    : invFactoryF_(invFA), invFactoryS_(invS), massInverseType_(Diagonal), schurCompOrdering_(false)
00092 {
00093    pcdParams_ = Teuchos::rcp(new Teuchos::ParameterList);
00094    lapParams_ = Teuchos::rcp(new Teuchos::ParameterList);
00095 
00096    lapParams_->set("Name",getPressureLaplaceString());
00097    pcdParams_->set("Name",getPCDString());
00098 
00099    buildTimers();
00100 }
00101 
00103 const Teko::LinearOp
00104 PCDStrategy::getHatInvA00(const Teko::BlockedLinearOp & A,BlockPreconditionerState & state) const
00105 {
00106    initializeState(A,state);
00107 
00108    return state.getModifiableOp("invF");
00109 }
00110 
00112 const Teko::LinearOp
00113 PCDStrategy::getTildeInvA00(const Teko::BlockedLinearOp & A,BlockPreconditionerState & state) const
00114 {
00115    initializeState(A,state);
00116 
00117    return state.getModifiableOp("invF");
00118 }
00119 
00121 const Teko::LinearOp
00122 PCDStrategy::getInvS(const Teko::BlockedLinearOp & A,BlockPreconditionerState & state) const
00123 {
00124    initializeState(A,state);
00125 
00126    return state.getLinearOp("invS");
00127 }
00128 
00129 void PCDStrategy::initializeState(const Teko::BlockedLinearOp & A,BlockPreconditionerState & state) const
00130 {
00131    Teko_DEBUG_SCOPE("PCDStrategy::initializeState",10);
00132    TEUCHOS_ASSERT(getRequestHandler()!=Teuchos::null);
00133 
00134    std::string pcdStr      = getPCDString();
00135    std::string presLapStr  = getPressureLaplaceString();
00136    std::string presMassStr = getPressureMassString();
00137 
00138    // no work to be done
00139    if(state.isInitialized())
00140       return;
00141 
00142    Teuchos::TimeMonitor timer(*initTimer_,true);
00143 
00144    // extract sub blocks
00145    LinearOp F  = Teko::getBlock(0,0,A);
00146    LinearOp Bt = Teko::getBlock(0,1,A);
00147    LinearOp B  = Teko::getBlock(1,0,A);
00148    LinearOp C  = Teko::getBlock(1,1,A);
00149 
00150    LinearOp  Qp = getRequestHandler()->request<LinearOp>(presMassStr);
00151    TEUCHOS_ASSERT(Qp!=Teuchos::null);
00152 
00153    // build the inverse Laplacian complement
00155    LinearOp iQp;
00156    if(massInverseType_==NotDiag) {
00157       ModifiableLinearOp & invMass = state.getModifiableOp("invMass");
00158       Teko_DEBUG_SCOPE("Building inv(Mass)",10);
00159 
00160       if(invMass==Teuchos::null)
00161          invMass = buildInverse(*invFactoryS_,Qp);
00162       else
00163          rebuildInverse(*invFactoryS_,Qp,invMass);
00164 
00165       iQp = invMass;
00166    }
00167    else {
00168       Teko_DEBUG_MSG("Building inverse mass of type \"" << Teko::getDiagonalName(massInverseType_) << "\"",10);
00169       iQp = getInvDiagonalOp(Qp,massInverseType_);
00170    }
00171 
00172    // build the inverse Laplacian complement
00174    ModifiableLinearOp & invLaplace = state.getModifiableOp("invLaplace");
00175    {
00176       Teuchos::TimeMonitor timer(*invSTimer_,true);
00177 
00178       // LinearOp laplace = getRequestHandler()->request<Teko::LinearOp>(presLapStr);
00179       LinearOp laplace = getRequestHandler()->request<Teko::LinearOp>(RequestMesg(lapParams_));
00180       TEUCHOS_ASSERT(laplace!=Teuchos::null);
00181       if(invLaplace==Teuchos::null)
00182          invLaplace = buildInverse(*invFactoryS_,laplace);
00183       else
00184          rebuildInverse(*invFactoryS_,laplace,invLaplace);
00185    }
00186 
00187    // build the inverse Schur complement
00189    {
00190       Teko_DEBUG_SCOPE("Building S",10);
00191       Teuchos::TimeMonitor timer(*opsTimer_,true);
00192 
00193       // build Schur-complement
00194       // LinearOp pcd = getRequestHandler()->request<Teko::LinearOp>(pcdStr);
00195       LinearOp pcd = getRequestHandler()->request<Teko::LinearOp>(RequestMesg(pcdParams_));
00196       TEUCHOS_ASSERT(pcd!=Teuchos::null);
00197       LinearOp invL = invLaplace;
00198  
00199       LinearOp invS;
00200       if(schurCompOrdering_==false)
00201          invS = multiply(iQp,pcd,invL);
00202       else
00203          invS = multiply(invL,pcd,iQp);
00204 
00205       state.addLinearOp("invS",invS);
00206    }
00207 
00208    // build inverse F
00210    {
00211       Teko_DEBUG_SCOPE("Building inv(F)",10);
00212       Teuchos::TimeMonitor timer(*invFTimer_,true);
00213 
00214       ModifiableLinearOp & invF = state.getModifiableOp("invF"); 
00215       if(invF==Teuchos::null)
00216          invF = buildInverse(*invFactoryF_,F);
00217       else
00218          rebuildInverse(*invFactoryF_,F,invF);
00219    }
00220 
00221    // mark state as initialized
00222    state.setInitialized(true);
00223 }
00224 
00236 void PCDStrategy::initializeFromParameterList(const Teuchos::ParameterList & pl,
00237                                               const InverseLibrary & invLib)
00238 {
00239    Teko_DEBUG_SCOPE("PCDStrategy::initializeFromParameterList",10);
00240 
00241    std::string invStr="Amesos", invFStr="", invSStr="";
00242    massInverseType_ = Diagonal;
00243 
00244    // "parse" the parameter list
00245    if(pl.isParameter("Inverse Type"))
00246       invStr = pl.get<std::string>("Inverse Type");
00247    if(pl.isParameter("Inverse F Type"))
00248       invFStr = pl.get<std::string>("Inverse F Type");
00249    if(pl.isParameter("Inverse Laplace Type"))
00250       invSStr = pl.get<std::string>("Inverse Laplace Type");
00251    if(pl.isParameter("Inverse Mass Type")) {
00252       std::string massInverseStr = pl.get<std::string>("Inverse Mass Type");
00253 
00254       // build inverse types
00255       massInverseType_ = getDiagonalType(massInverseStr);
00256    }
00257    if(pl.isParameter("Flip Schur Complement Ordering"))
00258       schurCompOrdering_ = pl.get<bool>("Flip Schur Complement Ordering");
00259 
00260    // set defaults as needed
00261    if(invFStr=="") invFStr = invStr;
00262    if(invSStr=="") invSStr = invStr;
00263 
00264    // read pressure laplace parameters
00265    if(pl.isSublist("Pressure Laplace Parameters"))
00266       lapParams_ = Teuchos::rcp(new Teuchos::ParameterList(pl.sublist("Pressure Laplace Parameters")));
00267    else
00268       lapParams_ = Teuchos::rcp(new Teuchos::ParameterList);
00269 
00270    // read pcd operator parameters
00271    if(pl.isSublist("Pressure Convection Diffusion Parameters"))
00272       pcdParams_ = Teuchos::rcp(new Teuchos::ParameterList(pl.sublist("Pressure Convection Diffusion Parameters")));
00273    else
00274       pcdParams_ = Teuchos::rcp(new Teuchos::ParameterList);
00275 
00276    // The user should not have already added this parameters
00277    TEUCHOS_TEST_FOR_EXCEPTION(lapParams_->isParameter("Name"),std::logic_error,
00278                    "Teko: Parameter \"Name\" is not allowed in the sublist \""+lapParams_->name()+"\"");
00279    TEUCHOS_TEST_FOR_EXCEPTION(lapParams_->isParameter("Tag"),std::logic_error,
00280                    "Teko: Parameter \"Tag\" is not allowed in the sublist \""+lapParams_->name()+"\"");
00281    TEUCHOS_TEST_FOR_EXCEPTION(pcdParams_->isParameter("Name"),std::logic_error,
00282                    "Teko: Parameter \"Name\" is not allowed in the sublist \""+pcdParams_->name()+"\"");
00283    TEUCHOS_TEST_FOR_EXCEPTION(pcdParams_->isParameter("Tag"),std::logic_error,
00284                    "Teko: Parameter \"Tag\" is not allowed in the sublist \""+pcdParams_->name()+"\"");
00285 
00286    Teko_DEBUG_MSG_BEGIN(5)
00287       DEBUG_STREAM << "PCD Strategy Parameters: " << std::endl;
00288       DEBUG_STREAM << "   inv type   = \"" << invStr  << "\"" << std::endl;
00289       DEBUG_STREAM << "   inv F type = \"" << invFStr << "\"" << std::endl;
00290       DEBUG_STREAM << "   inv Laplace type = \"" << invSStr << "\"" << std::endl;
00291       DEBUG_STREAM << "   inv Mass type = \"" << Teko::getDiagonalName(massInverseType_) << "\"" << std::endl;
00292       DEBUG_STREAM << "PCD Strategy Parameter list: " << std::endl;
00293       pl.print(DEBUG_STREAM);
00294    Teko_DEBUG_MSG_END()
00295 
00296    // build velocity inverse factory
00297    invFactoryF_ = invLib.getInverseFactory(invFStr);
00298  
00299    if(invFStr==invSStr)
00300       invFactoryS_ = invFactoryF_;
00301    else
00302       invFactoryS_ = invLib.getInverseFactory(invSStr);
00303 
00304    lapParams_->set("Name",getPressureLaplaceString());
00305    pcdParams_->set("Name",getPCDString());
00306 
00307    // setup a request for required operators
00308    getRequestHandler()->preRequest<Teko::LinearOp>(getPressureMassString());
00309    // getRequestHandler()->preRequest<Teko::LinearOp>(getPCDString());
00310    // getRequestHandler()->preRequest<Teko::LinearOp>(getPressureLaplaceString());
00311    getRequestHandler()->preRequest<Teko::LinearOp>(Teko::RequestMesg(lapParams_));
00312    getRequestHandler()->preRequest<Teko::LinearOp>(Teko::RequestMesg(pcdParams_));
00313 }
00314 
00316 Teuchos::RCP<Teuchos::ParameterList> PCDStrategy::getRequestedParameters() const 
00317 {
00318    TEUCHOS_ASSERT(false);
00319 
00320    return Teuchos::null;
00321 }
00322 
00324 bool PCDStrategy::updateRequestedParameters(const Teuchos::ParameterList & pl) 
00325 {
00326    TEUCHOS_ASSERT(false);
00327 
00328    return true;
00329 }
00330 
00331 } // end namespace NS
00332 } // end namespace Teko
 All Classes Files Functions Variables