Teko Version of the Day
Teko_LSCPreconditionerFactory.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_LSCPreconditionerFactory.hpp"
00048 
00049 #include "Thyra_DefaultMultipliedLinearOp.hpp"
00050 #include "Thyra_DefaultAddedLinearOp.hpp"
00051 #include "Thyra_DefaultIdentityLinearOp.hpp"
00052 #include "Thyra_DefaultZeroLinearOp.hpp"
00053 #include "Thyra_get_Epetra_Operator.hpp"
00054 
00055 #include "Teko_LU2x2InverseOp.hpp"
00056 #include "Teko_Utilities.hpp"
00057 #include "Teko_BlockUpperTriInverseOp.hpp"
00058 #include "Teko_StaticLSCStrategy.hpp"
00059 #include "Teko_InvLSCStrategy.hpp"
00060 #include "Teko_PresLaplaceLSCStrategy.hpp"
00061 #include "Teko_LSCSIMPLECStrategy.hpp"
00062 
00063 #include "EpetraExt_RowMatrixOut.h"
00064 
00065 #include "Teuchos_Time.hpp"
00066 
00067 namespace Teko {
00068 namespace NS {
00069 
00070 using Teuchos::rcp;
00071 using Teuchos::rcp_dynamic_cast;
00072 using Teuchos::RCP;
00073 
00074 using Thyra::multiply;
00075 using Thyra::add;
00076 using Thyra::identity;
00077 
00078 // Stabilized constructor
00079 LSCPreconditionerFactory::LSCPreconditionerFactory(const LinearOp & invF,const LinearOp & invBQBtmC,
00080                                                    const LinearOp & invD,const LinearOp & invMass)
00081       : invOpsStrategy_(rcp(new StaticLSCStrategy(invF,invBQBtmC,invD,invMass))), isSymmetric_(true)
00082 { }
00083 
00084 // Stable constructor
00085 LSCPreconditionerFactory::LSCPreconditionerFactory(const LinearOp & invF, const LinearOp & invBQBtmC,
00086                                                    const LinearOp & invMass)
00087       : invOpsStrategy_(rcp(new StaticLSCStrategy(invF,invBQBtmC,invMass))), isSymmetric_(true)
00088 { }
00089 
00090 // fully generic constructor
00091 LSCPreconditionerFactory::LSCPreconditionerFactory(const RCP<LSCStrategy> & strategy)
00092    : invOpsStrategy_(strategy), isSymmetric_(true)
00093 { }
00094 
00095 LSCPreconditionerFactory::LSCPreconditionerFactory() : isSymmetric_(true)
00096 { }
00097 
00098 // for PreconditionerFactoryBase
00100 
00101 // initialize a newly created preconditioner object
00102 LinearOp LSCPreconditionerFactory::buildPreconditionerOperator(BlockedLinearOp & blockOp,BlockPreconditionerState & state) const
00103 {
00104    Teko_DEBUG_SCOPE("LSCPreconditionerFactory::buildPreconditionerOperator",10);
00105    Teko_DEBUG_EXPR(Teuchos::Time timer(""));
00106    Teko_DEBUG_EXPR(Teuchos::Time totalTimer(""));
00107    Teko_DEBUG_EXPR(totalTimer.start());
00108 
00109    // extract sub-matrices from source operator 
00110    LinearOp F  = blockOp->getBlock(0,0);
00111    LinearOp B  = blockOp->getBlock(1,0);
00112    LinearOp Bt = blockOp->getBlock(0,1);
00113 
00114    if(not isSymmetric_)
00115       Bt = scale(-1.0,adjoint(B));
00116 
00117    // build what is neccessary for the state object
00118    Teko_DEBUG_EXPR(timer.start(true));
00119    invOpsStrategy_->buildState(blockOp,state);
00120    Teko_DEBUG_EXPR(timer.stop());
00121    Teko_DEBUG_MSG("LSCPrecFact::buildPO BuildStateTime = " << timer.totalElapsedTime(),2);
00122 
00123    // extract operators from strategy
00124    Teko_DEBUG_EXPR(timer.start(true));
00125    LinearOp invF      = invOpsStrategy_->getInvF(blockOp,state);
00126    LinearOp invBQBtmC = invOpsStrategy_->getInvBQBt(blockOp,state);
00127    LinearOp invBHBtmC = invOpsStrategy_->getInvBHBt(blockOp,state);
00128    LinearOp outerStab = invOpsStrategy_->getOuterStabilization(blockOp,state);
00129    LinearOp innerStab = invOpsStrategy_->getInnerStabilization(blockOp,state);
00130 
00131    // if necessary build an identity mass matrix
00132    LinearOp invMass   = invOpsStrategy_->getInvMass(blockOp,state);
00133    LinearOp HScaling  = invOpsStrategy_->getHScaling(blockOp,state);
00134    if(invMass==Teuchos::null)  invMass = identity<double>(F->range());
00135    if(HScaling==Teuchos::null) HScaling = identity<double>(F->range());
00136    Teko_DEBUG_EXPR(timer.stop());
00137    Teko_DEBUG_MSG("LSCPrecFact::buildPO GetInvTime = " << timer.totalElapsedTime(),2);
00138 
00139    // need to build Schur complement,  inv(P) = inv(B*Bt)*(B*F*Bt)*inv(B*Bt)
00140 
00141    // first construct middle operator: M = B * inv(Mass) * F * inv(Mass) * Bt
00142    LinearOp M = 
00143       //          (B * inv(Mass) ) * F * (inv(Mass) * Bt)
00144       multiply( multiply(B,invMass), F , multiply(HScaling,Bt));
00145    if(innerStab!=Teuchos::null) // deal with inner stabilization
00146       M = add(M, innerStab);
00147       
00148    // now construct a linear operator schur complement
00149    LinearOp invPschur; 
00150    if(outerStab!=Teuchos::null)
00151       invPschur = add(multiply(invBQBtmC, M , invBHBtmC), outerStab);
00152    else
00153       invPschur = multiply(invBQBtmC, M , invBHBtmC);
00154 
00155    // build the preconditioner operator: Use LDU or upper triangular approximation
00156    if(invOpsStrategy_->useFullLDU()) { 
00157       Teko_DEBUG_EXPR(totalTimer.stop());
00158       Teko_DEBUG_MSG("LSCPrecFact::buildPO TotalTime = " << totalTimer.totalElapsedTime(),2);
00159 
00160       // solve using a full LDU decomposition
00161       return createLU2x2InverseOp(blockOp,invF,invPschur,"LSC-LDU");
00162    } else {
00163       // build diagonal operations
00164       std::vector<LinearOp> invDiag(2);
00165       invDiag[0] = invF;
00166       invDiag[1] = Thyra::scale(-1.0,invPschur);
00167 
00168       // get upper triangular matrix
00169       BlockedLinearOp U = getUpperTriBlocks(blockOp); 
00170 
00171       Teko_DEBUG_EXPR(totalTimer.stop());
00172       Teko_DEBUG_MSG("LSCPrecFact::buildPO TotalTime = " << totalTimer.totalElapsedTime(),2);
00173 
00174       // solve using only one inversion of F
00175       return createBlockUpperTriInverseOp(U,invDiag,"LSC-Upper");
00176    }
00177 }
00178 
00180 void LSCPreconditionerFactory::initializeFromParameterList(const Teuchos::ParameterList & pl)
00181 {
00182    Teko_DEBUG_SCOPE("LSCPreconditionerFactory::initializeFromParameterList",10);
00183 
00184    RCP<const InverseLibrary> invLib = getInverseLibrary();
00185 
00186    if(pl.isParameter("Is Symmetric"))
00187       isSymmetric_ = pl.get<bool>("Is Symmetric");
00188 
00189    std::string name = "Basic Inverse";
00190    if(pl.isParameter("Strategy Name"))
00191       name = pl.get<std::string>("Strategy Name");
00192    const Teuchos::ParameterEntry * pe = pl.getEntryPtr("Strategy Settings");
00193 
00194    // check for a mistake in input file
00195    if(name!="Basic Inverse" && pe==0) {
00196       RCP<Teuchos::FancyOStream> out = getOutputStream();
00197       *out << "LSC Construction failed: ";
00198       *out << "Strategy \"" << name << "\" requires a \"Strategy Settings\" sublist" << std::endl;
00199       throw std::runtime_error("LSC Construction failed: Strategy Settings not set");
00200    }
00201 
00202    // get the parameter list to construct the strategy
00203    Teuchos::RCP<const Teuchos::ParameterList> stratPL = Teuchos::rcpFromRef(pl);
00204    if(pe!=0)
00205       stratPL = Teuchos::rcpFromRef(pl.sublist("Strategy Settings"));
00206 
00207    // build the strategy object
00208    RCP<LSCStrategy> strategy = buildStrategy(name,*stratPL,invLib,getRequestHandler());
00209  
00210    // strategy could not be built
00211    if(strategy==Teuchos::null) {
00212       RCP<Teuchos::FancyOStream> out = getOutputStream();
00213       *out << "LSC Construction failed: ";
00214       *out << "Strategy \"" << name << "\" could not be constructed" << std::endl;
00215       throw std::runtime_error("LSC Construction failed: Strategy could not be constructed");
00216    }
00217 
00218    strategy->setSymmetric(isSymmetric_);
00219    invOpsStrategy_ = strategy;
00220 }
00221 
00223 Teuchos::RCP<Teuchos::ParameterList> LSCPreconditionerFactory::getRequestedParameters() const
00224 {
00225    return invOpsStrategy_->getRequestedParameters();
00226 }
00227 
00229 bool LSCPreconditionerFactory::updateRequestedParameters(const Teuchos::ParameterList & pl)
00230 {
00231    return invOpsStrategy_->updateRequestedParameters(pl);
00232 }
00233 
00235 // Static members and methods
00237 
00239 CloneFactory<LSCStrategy> LSCPreconditionerFactory::strategyBuilder_;
00240 
00253 RCP<LSCStrategy> LSCPreconditionerFactory::buildStrategy(const std::string & name, 
00254                                                          const Teuchos::ParameterList & settings,
00255                                                          const RCP<const InverseLibrary> & invLib,
00256                                                          const RCP<RequestHandler> & rh)
00257 {
00258    Teko_DEBUG_SCOPE("LSCPreconditionerFactory::buildStrategy",10);
00259 
00260    // initialize the defaults if necessary
00261    if(strategyBuilder_.cloneCount()==0) initializeStrategyBuilder();
00262 
00263    // request the preconditioner factory from the CloneFactory
00264    Teko_DEBUG_MSG("Building LSC strategy \"" << name << "\"",1);
00265    RCP<LSCStrategy> strategy = strategyBuilder_.build(name);
00266 
00267    if(strategy==Teuchos::null) return Teuchos::null;
00268 
00269    // now that inverse library has been set,
00270    // pass in the parameter list
00271    strategy->setRequestHandler(rh);
00272    strategy->initializeFromParameterList(settings,*invLib);
00273 
00274    return strategy;
00275 }
00276 
00290 void LSCPreconditionerFactory::addStrategy(const std::string & name,const RCP<Cloneable> & clone)
00291 {
00292    // initialize the defaults if necessary
00293    if(strategyBuilder_.cloneCount()==0) initializeStrategyBuilder();
00294 
00295    // add clone to builder
00296    strategyBuilder_.addClone(name,clone); 
00297 }
00298 
00300 void LSCPreconditionerFactory::initializeStrategyBuilder()
00301 {
00302    RCP<Cloneable> clone;
00303 
00304    // add various strategies to the factory
00305    clone = rcp(new AutoClone<InvLSCStrategy>());
00306    strategyBuilder_.addClone("Basic Inverse",clone);
00307 
00308    // add various strategies to the factory
00309    clone = rcp(new AutoClone<PresLaplaceLSCStrategy>());
00310    strategyBuilder_.addClone("Pressure Laplace",clone);
00311 
00312    // add various strategies to the factory
00313    clone = rcp(new AutoClone<LSCSIMPLECStrategy>());
00314    strategyBuilder_.addClone("SIMPLEC",clone);
00315 }
00316 
00317 } // end namespace NS
00318 } // end namespace Teko
 All Classes Files Functions Variables