Teko Version of the Day
Teko_ModALPreconditionerFactory.cpp
00001 /*
00002  * Author: Zhen Wang
00003  * Email: wangz@ornl.gov
00004  *        zhen.wang@alum.emory.edu
00005  */
00006 
00007 #include "Teko_ModALPreconditionerFactory.hpp"
00008 
00009 #include "Thyra_DefaultMultipliedLinearOp.hpp"
00010 #include "Thyra_DefaultAddedLinearOp.hpp"
00011 #include "Thyra_DefaultIdentityLinearOp.hpp"
00012 #include "Thyra_DefaultZeroLinearOp.hpp"
00013 #include "Thyra_get_Epetra_Operator.hpp"
00014 
00015 #include "Teko_LU2x2InverseOp.hpp"
00016 #include "Teko_Utilities.hpp"
00017 #include "Teko_BlockLowerTriInverseOp.hpp"
00018 #include "Teko_BlockUpperTriInverseOp.hpp"
00019 #include "Teko_StaticLSCStrategy.hpp"
00020 #include "Teko_InvLSCStrategy.hpp"
00021 #include "Teko_PresLaplaceLSCStrategy.hpp"
00022 
00023 #include "EpetraExt_RowMatrixOut.h"
00024 
00025 #include "Teuchos_Time.hpp"
00026 
00027 namespace Teko
00028 {
00029 
00030 namespace NS
00031 {
00032 
00033 ModALPrecondState::ModALPrecondState():
00034 pressureMassMatrix_(Teuchos::null), invPressureMassMatrix_(Teuchos::null)
00035 {
00036 }
00037 
00038 ModALPreconditionerFactory::ModALPreconditionerFactory(const Teuchos::RCP<InverseFactory> & factory) :
00039       invOpsStrategy_(Teuchos::rcp(new InvModALStrategy(factory))),
00040       isSymmetric_(true)
00041 {
00042 }
00043 
00044 ModALPreconditionerFactory::ModALPreconditionerFactory(const Teuchos::RCP<InverseFactory> & invFactoryA,
00045       const Teuchos::RCP<InverseFactory> & invFactoryS) :
00046       invOpsStrategy_(Teuchos::rcp(new InvModALStrategy(invFactoryA, invFactoryS))),
00047       isSymmetric_(true)
00048 {
00049 }
00050 
00051 ModALPreconditionerFactory::ModALPreconditionerFactory(const Teuchos::RCP<InverseFactory> & factory,
00052       LinearOp & pressureMassMatrix) :
00053    invOpsStrategy_(Teuchos::rcp(new InvModALStrategy(factory, pressureMassMatrix))), isSymmetric_(true)
00054 {
00055 }
00056 
00057 ModALPreconditionerFactory::ModALPreconditionerFactory(const Teuchos::RCP<InverseFactory> & invFactoryA,
00058       const Teuchos::RCP<InverseFactory> & invFactoryS,
00059       LinearOp & pressureMassMatrix) :
00060    invOpsStrategy_(Teuchos::rcp(new InvModALStrategy(invFactoryA, invFactoryS, pressureMassMatrix))),
00061    isSymmetric_(true)
00062 {
00063 }
00064 
00065 // Construct from a strategy.
00066 ModALPreconditionerFactory::ModALPreconditionerFactory(const Teuchos::RCP<InvModALStrategy> & strategy) :
00067    invOpsStrategy_(strategy), isSymmetric_(true)
00068 {
00069 }
00070 
00071 LinearOp ModALPreconditionerFactory::buildPreconditionerOperator(BlockedLinearOp & alOp,
00072       BlockPreconditionerState & state) const
00073 {
00074    Teko_DEBUG_SCOPE("ModALPreconditionerFactory::buildPreconditionerOperator()", 10);
00075    Teko_DEBUG_EXPR(Teuchos::Time timer(""));
00076    Teko_DEBUG_EXPR(Teuchos::Time totalTimer(""));
00077    Teko_DEBUG_EXPR(totalTimer.start());
00078 
00079    // Only for 2D or 3D problems.
00080    int dim = blockRowCount(alOp) - 1;
00081    TEUCHOS_ASSERT(dim == 2 || dim == 3);
00082 
00083    // Build what is necessary for the state object.
00084    Teko_DEBUG_EXPR(timer.start(true));
00085    invOpsStrategy_->buildState(alOp, state);
00086    Teko_DEBUG_EXPR(timer.stop());
00087    Teko_DEBUG_MSG("ModALPreconditionerFactory::buildPreconditionerOperator():BuildStateTime = "
00088          << timer.totalElapsedTime(), 2);
00089 
00090    // Extract inverse operators from strategy
00091    Teko_DEBUG_EXPR(timer.start(true));
00092    LinearOp invA11p = invOpsStrategy_->getInvA11p(state);
00093    LinearOp invA22p = invOpsStrategy_->getInvA22p(state);
00094    LinearOp invA33p;
00095    if(dim == 3)
00096    {
00097       invA33p = invOpsStrategy_->getInvA33p(state);
00098    }
00099 
00100    // The inverse of S can be built from strategy,
00101    // or just a diagonal matrix.
00102    ModALPrecondState * modALState = dynamic_cast<ModALPrecondState*>(&state);
00103    TEUCHOS_ASSERT(modALState != NULL);
00104    LinearOp invS;
00105    if(modALState->isStabilized_)
00106    {
00107       invS = invOpsStrategy_->getInvS(state);
00108    }
00109    else
00110    {
00111       invS = scale(modALState->gamma_, modALState->invPressureMassMatrix_);
00112    }
00113 
00114    Teko_DEBUG_EXPR(timer.stop());
00115    Teko_DEBUG_MSG("ModALPrecFact::buildPreconditionerOperator(): GetInvTime = "
00116          << timer.totalElapsedTime(), 2);
00117 
00118    // Build diagonal operations.
00119    std::vector<LinearOp> invDiag;
00120    invDiag.resize(dim + 1);
00121    invDiag[0] = invA11p;
00122    invDiag[1] = invA22p;
00123    if(dim == 2)
00124    {
00125       invDiag[2] = scale(-1.0, invS);
00126    }
00127    else if(dim == 3)
00128    {
00129       invDiag[2] = invA33p;
00130       invDiag[3] = scale(-1.0, invS);
00131    }
00132 
00133    // Get the upper triangular matrix.
00134    BlockedLinearOp U = getUpperTriBlocks(alOp);
00135 
00136    Teko_DEBUG_EXPR(totalTimer.stop());
00137    Teko_DEBUG_MSG("ModALPrecFact::buildPreconditionerOperator TotalTime = "
00138          << totalTimer.totalElapsedTime(), 2);
00139 
00140    // Create the preconditioner
00141    return createBlockUpperTriInverseOp(U, invDiag, "Modified AL preconditioner-Upper");
00142 }
00143 
00144 } // end namespace NS
00145 
00146 } // end namespace Teko
 All Classes Files Functions Variables