Teko Version of the Day
Teko_LSCSIMPLECStrategy.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_LSCSIMPLECStrategy.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 "Teko_LSCPreconditionerFactory.hpp"
00068 #include "Teko_EpetraHelpers.hpp"
00069 #include "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 // LSCSIMPLECStrategy Implementation
00081 
00082 // constructors
00084 LSCSIMPLECStrategy::LSCSIMPLECStrategy()
00085    : invFactoryF_(Teuchos::null), invFactoryS_(Teuchos::null)
00086    , useFullLDU_(false), scaleType_(Diagonal)
00087 { }
00088 
00090 
00091 void LSCSIMPLECStrategy::buildState(BlockedLinearOp & A,BlockPreconditionerState & state) const
00092 {
00093    Teko_DEBUG_SCOPE("LSCSIMPLECStrategy::buildState",10);
00094 
00095    LSCPrecondState * lscState = dynamic_cast<LSCPrecondState*>(&state);
00096    TEUCHOS_ASSERT(lscState!=0);
00097 
00098    // if neccessary save state information
00099    if(not lscState->isInitialized()) {
00100       Teko_DEBUG_EXPR(Teuchos::Time timer(""));
00101 
00102       // construct operators
00103       {
00104          Teko_DEBUG_SCOPE("LSC-SIMPLEC::buildState constructing operators",1);
00105          Teko_DEBUG_EXPR(timer.start(true));
00106 
00107          initializeState(A,lscState);
00108 
00109          Teko_DEBUG_EXPR(timer.stop());
00110          Teko_DEBUG_MSG("LSC-SIMPLEC::buildState BuildOpsTime = " << timer.totalElapsedTime(),1);
00111       }
00112 
00113       // Build the inverses
00114       {
00115          Teko_DEBUG_SCOPE("LSC-SIMPLEC::buildState calculating inverses",1);
00116          Teko_DEBUG_EXPR(timer.start(true));
00117 
00118          computeInverses(A,lscState);
00119 
00120          Teko_DEBUG_EXPR(timer.stop());
00121          Teko_DEBUG_MSG("LSC-SIMPLEC::buildState BuildInvTime = " << timer.totalElapsedTime(),1);
00122       }
00123    }
00124 }
00125 
00126 // functions inherited from LSCStrategy
00127 LinearOp LSCSIMPLECStrategy::getInvBQBt(const BlockedLinearOp & A,BlockPreconditionerState & state) const
00128 {
00129    return state.getInverse("invBQBtmC");
00130 }
00131 
00132 LinearOp LSCSIMPLECStrategy::getInvBHBt(const BlockedLinearOp & A,BlockPreconditionerState & state) const
00133 {
00134    return state.getInverse("invBQBtmC").getConst();
00135 }
00136 
00137 LinearOp LSCSIMPLECStrategy::getInvF(const BlockedLinearOp & A,BlockPreconditionerState & state) const
00138 {
00139    return state.getInverse("invF");
00140 }
00141 
00142 LinearOp LSCSIMPLECStrategy::getInnerStabilization(const BlockedLinearOp & A,BlockPreconditionerState & state) const
00143 {
00144    LSCPrecondState * lscState = dynamic_cast<LSCPrecondState*>(&state);
00145    TEUCHOS_ASSERT(lscState!=0);
00146    TEUCHOS_ASSERT(lscState->isInitialized())
00147 
00148    const LinearOp C  = getBlock(1,1,A);
00149    return scale(-1.0,C);
00150 }
00151 
00152 
00153 LinearOp LSCSIMPLECStrategy::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 LSCSIMPLECStrategy::getHScaling(const BlockedLinearOp & A,BlockPreconditionerState & state) const
00163 {
00164    return getInvMass(A,state);
00165 }
00166 
00168 void LSCSIMPLECStrategy::initializeState(const BlockedLinearOp & A,LSCPrecondState * state) const
00169 {
00170    Teko_DEBUG_SCOPE("LSCSIMPLECStrategy::initializeState",10);
00171 
00172    const LinearOp F  = getBlock(0,0,A);
00173    const LinearOp Bt = getBlock(0,1,A);
00174    const LinearOp B  = getBlock(1,0,A);
00175    const LinearOp C  = getBlock(1,1,A);
00176 
00177    bool isStabilized = (not isZeroOp(C));
00178 
00179    state->invMass_ = getInvDiagonalOp(F,scaleType_);
00180 
00181    // compute BQBt
00182    state->BQBt_ = explicitMultiply(B,state->invMass_,Bt,state->BQBt_);
00183    if(isStabilized) {
00184       // now build B*Q*Bt-C
00185       Teko::ModifiableLinearOp BQBtmC = state->getInverse("BQBtmC");
00186       BQBtmC = explicitAdd(state->BQBt_,scale(-1.0,C),BQBtmC);
00187       state->addInverse("BQBtmC",BQBtmC);
00188    }
00189    Teko_DEBUG_MSG("Computed BQBt",10);
00190 
00191    state->setInitialized(true);
00192 }
00193 
00199 void LSCSIMPLECStrategy::computeInverses(const BlockedLinearOp & A,LSCPrecondState * state) const
00200 {
00201    Teko_DEBUG_SCOPE("LSCSIMPLECStrategy::computeInverses",10);
00202    Teko_DEBUG_EXPR(Teuchos::Time invTimer(""));
00203 
00204    const LinearOp F  = getBlock(0,0,A);
00205 
00207 
00208    // (re)build the inverse of F
00209    Teko_DEBUG_MSG("LSC-SIMPLEC::computeInverses Building inv(F)",1);
00210    Teko_DEBUG_EXPR(invTimer.start(true));
00211    InverseLinearOp invF = state->getInverse("invF");
00212    if(invF==Teuchos::null) {
00213       invF = buildInverse(*invFactoryF_,F);
00214       state->addInverse("invF",invF); 
00215    } else {
00216       rebuildInverse(*invFactoryF_,F,invF);
00217    }
00218    Teko_DEBUG_EXPR(invTimer.stop());
00219    Teko_DEBUG_MSG("LSC-SIMPLEC::computeInverses GetInvF = " << invTimer.totalElapsedTime(),1);
00220 
00222 
00223    // (re)build the inverse of BQBt 
00224    Teko_DEBUG_MSG("LSC-SIMPLEC::computeInverses Building inv(BQBtmC)",1);
00225    Teko_DEBUG_EXPR(invTimer.start(true));
00226    const LinearOp BQBt = state->getInverse("BQBtmC");
00227    InverseLinearOp invBQBt = state->getInverse("invBQBtmC");
00228    if(invBQBt==Teuchos::null) {
00229       invBQBt = buildInverse(*invFactoryS_,BQBt);
00230       state->addInverse("invBQBtmC",invBQBt); 
00231    } else {
00232       rebuildInverse(*invFactoryS_,BQBt,invBQBt);
00233    }
00234    Teko_DEBUG_EXPR(invTimer.stop());
00235    Teko_DEBUG_MSG("LSC-SIMPLEC::computeInverses GetInvBQBt = " << invTimer.totalElapsedTime(),1);
00236 }
00237 
00239 void LSCSIMPLECStrategy::initializeFromParameterList(const Teuchos::ParameterList & pl,const InverseLibrary & invLib) 
00240 {
00241    // get string specifying inverse
00242    std::string invStr="", invVStr="", invPStr="";
00243    bool useLDU = false;
00244    scaleType_ = Diagonal;
00245 
00246    // "parse" the parameter list
00247    if(pl.isParameter("Inverse Type"))
00248       invStr = pl.get<std::string>("Inverse Type");
00249    if(pl.isParameter("Inverse Velocity Type"))
00250       invVStr = pl.get<std::string>("Inverse Velocity Type");
00251    if(pl.isParameter("Inverse Pressure Type")) 
00252       invPStr = pl.get<std::string>("Inverse Pressure Type");
00253    if(pl.isParameter("Use LDU"))
00254       useLDU = pl.get<bool>("Use LDU");
00255    if(pl.isParameter("Scaling Type")) {
00256       scaleType_ = getDiagonalType(pl.get<std::string>("Scaling Type"));
00257       TEUCHOS_TEST_FOR_EXCEPT(scaleType_==NotDiag);
00258    }
00259 
00260    Teko_DEBUG_MSG_BEGIN(0)
00261       DEBUG_STREAM << "LSC Inverse Strategy Parameters: " << std::endl;
00262       DEBUG_STREAM << "   inv type   = \"" << invStr  << "\"" << std::endl;
00263       DEBUG_STREAM << "   inv v type = \"" << invVStr << "\"" << std::endl;
00264       DEBUG_STREAM << "   inv p type = \"" << invPStr << "\"" << std::endl;
00265       DEBUG_STREAM << "   use ldu    = " << useLDU << std::endl;
00266       DEBUG_STREAM << "   scale type    = " << getDiagonalName(scaleType_) << std::endl;
00267       DEBUG_STREAM << "LSC  Inverse Strategy Parameter list: " << std::endl;
00268       pl.print(DEBUG_STREAM);
00269    Teko_DEBUG_MSG_END()
00270 
00271    // set defaults as needed
00272    if(invStr=="") invStr = "Amesos";
00273    if(invVStr=="") invVStr = invStr;
00274    if(invPStr=="") invPStr = invStr;
00275 
00276    // build velocity inverse factory
00277    invFactoryF_ = invLib.getInverseFactory(invVStr);
00278    invFactoryS_ = invFactoryF_; // by default these are the same
00279    if(invVStr!=invPStr) // if different, build pressure inverse factory
00280       invFactoryS_ = invLib.getInverseFactory(invPStr);
00281 
00282    // set other parameters
00283    setUseFullLDU(useLDU);
00284 }
00285 
00286 } // end namespace NS
00287 } // end namespace Teko
 All Classes Files Functions Variables