Teko Version of the Day
Teko_InvModALStrategy.cpp
00001 /*
00002  * Author: Zhen Wang
00003  * Email: wangz@ornl.gov
00004  *        zhen.wang@alum.emory.edu
00005  */
00006 
00007 #include "Thyra_DefaultDiagonalLinearOp.hpp"
00008 #include "Thyra_EpetraThyraWrappers.hpp"
00009 #include "Thyra_get_Epetra_Operator.hpp"
00010 #include "Thyra_EpetraLinearOp.hpp"
00011 #include "Thyra_DefaultAddedLinearOp_def.hpp"
00012 #include "Thyra_DefaultIdentityLinearOp_decl.hpp"
00013 
00014 #include "Epetra_Vector.h"
00015 #include "Epetra_Map.h"
00016 
00017 #include "EpetraExt_RowMatrixOut.h"
00018 #include "EpetraExt_MultiVectorOut.h"
00019 
00020 #include "Teuchos_Time.hpp"
00021 
00022 #include "Teko_Utilities.hpp"
00023 
00024 #include "Teko_InvModALStrategy.hpp"
00025 #include "Teko_ModALPreconditionerFactory.hpp"
00026 
00027 using Teuchos::RCP;
00028 using Teuchos::rcp_dynamic_cast;
00029 using Teuchos::rcp_const_cast;
00030 
00031 namespace Teko
00032 {
00033 
00034 namespace NS
00035 {
00036 
00037 // Empty constructor.
00038 InvModALStrategy::InvModALStrategy() :
00039       invFactoryA_(Teuchos::null), invFactoryS_(Teuchos::null),
00040       pressureMassMatrix_(Teuchos::null), gamma_(0.05),
00041       scaleType_(Diagonal), isSymmetric_(true)
00042 {
00043 }
00044 
00045 // If there is only one InverseFactory, use it for all solves.
00046 InvModALStrategy::InvModALStrategy(const Teuchos::RCP<InverseFactory> & factory) :
00047       invFactoryA_(factory), invFactoryS_(factory),
00048       pressureMassMatrix_(Teuchos::null), gamma_(0.05),
00049       scaleType_(Diagonal), isSymmetric_(true)
00050 {
00051 }
00052 
00053 // If there are two InverseFactory's...
00054 InvModALStrategy::InvModALStrategy(const Teuchos::RCP<InverseFactory> & invFactoryA,
00055       const Teuchos::RCP<InverseFactory> & invFactoryS) :
00056       invFactoryA_(invFactoryA), invFactoryS_(invFactoryS),
00057       pressureMassMatrix_(Teuchos::null), gamma_(0.05),
00058       scaleType_(Diagonal), isSymmetric_(true)
00059 {
00060 }
00061 
00062 // If there are two InverseFactory's...
00063 InvModALStrategy::InvModALStrategy(const Teuchos::RCP<InverseFactory> & invFactory,
00064       LinearOp & pressureMassMatrix) :
00065       invFactoryA_(invFactory), invFactoryS_(invFactory),
00066       pressureMassMatrix_(pressureMassMatrix), gamma_(0.05),
00067       scaleType_(Diagonal), isSymmetric_(true)
00068 {
00069 }
00070 
00071 // If there are two InverseFactory's...
00072 InvModALStrategy::InvModALStrategy(const Teuchos::RCP<InverseFactory> & invFactoryA,
00073       const Teuchos::RCP<InverseFactory> & invFactoryS, LinearOp & pressureMassMatrix) :
00074       invFactoryA_(invFactoryA), invFactoryS_(invFactoryS),
00075       pressureMassMatrix_(pressureMassMatrix), gamma_(0.05),
00076       scaleType_(Diagonal), isSymmetric_(true)
00077 {
00078 }
00079 
00080 // Return "inverses".
00081 LinearOp InvModALStrategy::getInvA11p(BlockPreconditionerState & state) const
00082 {
00083    return state.getInverse("invA11p");
00084 }
00085 
00086 LinearOp InvModALStrategy::getInvA22p(BlockPreconditionerState & state) const
00087 {
00088    return state.getInverse("invA22p");
00089 }
00090 
00091 LinearOp InvModALStrategy::getInvA33p(BlockPreconditionerState & state) const
00092 {
00093    return state.getInverse("invA33p");
00094 }
00095 
00096 LinearOp InvModALStrategy::getInvS(BlockPreconditionerState & state) const
00097 {
00098    return state.getInverse("invS");
00099 }
00100 
00101 // Set pressure mass matrix.
00102 void InvModALStrategy::setPressureMassMatrix(const LinearOp & pressureMassMatrix)
00103 {
00104    pressureMassMatrix_ = pressureMassMatrix;
00105 }
00106 
00107 // Set gamma.
00108 void InvModALStrategy::setGamma(double gamma)
00109 {
00110    TEUCHOS_ASSERT(gamma > 0.0);
00111    gamma_ = gamma;
00112 }
00113 
00114 void InvModALStrategy::buildState(const BlockedLinearOp & alOp,
00115       BlockPreconditionerState & state) const
00116 {
00117    Teko_DEBUG_SCOPE("InvModALStrategy::buildState", 10);
00118 
00119    ModALPrecondState * modALState = dynamic_cast<ModALPrecondState*>(&state);
00120    TEUCHOS_ASSERT(modALState != NULL);
00121 
00122    // if necessary save state information
00123    if(not modALState->isInitialized())
00124    {
00125       Teko_DEBUG_EXPR(Teuchos::Time timer(""));
00126 
00127       {
00128          // construct operators
00129          Teko_DEBUG_SCOPE("ModAL::buildState: Initializing state object", 1);
00130          Teko_DEBUG_EXPR(timer.start(true));
00131 
00132          initializeState(alOp, modALState);
00133 
00134          Teko_DEBUG_EXPR(timer.stop());
00135          Teko_DEBUG_MSG("ModAL::buildState: BuildOpsTime = " << timer.totalElapsedTime(), 1);
00136       }
00137 
00138       {
00139          // Build the inverses
00140          Teko_DEBUG_SCOPE("ModAL::buildState: Computing inverses", 1);
00141          Teko_DEBUG_EXPR(timer.start(true));
00142 
00143          computeInverses(alOp, modALState);
00144 
00145          Teko_DEBUG_EXPR(timer.stop());
00146          Teko_DEBUG_MSG("ModAL::buildState: BuildInvTime = " << timer.totalElapsedTime(), 1);
00147       }
00148    }
00149 }
00150 
00151 // Initialize the state object using the ALOperator.
00152 void InvModALStrategy::initializeState(const BlockedLinearOp & alOp,
00153       ModALPrecondState *state) const
00154 {
00155    Teko_DEBUG_SCOPE("InvModALStrategy::initializeState", 10);
00156 
00157    // Extract sub-matrices from blocked linear operator.
00158    int dim = blockRowCount(alOp) - 1;
00159    TEUCHOS_ASSERT(dim == 2 || dim == 3);
00160 
00161    LinearOp lpA11 = getBlock(0, 0, alOp);
00162    LinearOp lpA22 = getBlock(1, 1, alOp);
00163    LinearOp lpA33, lpB1, lpB2, lpB3, lpB1t, lpB2t, lpB3t, lpC;
00164 
00165    // 2D problem.
00166    if(dim == 2)
00167    {
00168       lpB1  = getBlock(2, 0, alOp);
00169       lpB2  = getBlock(2, 1, alOp);
00170       lpB1t = getBlock(0, 2, alOp);
00171       lpB2t = getBlock(1, 2, alOp);
00172       lpC = getBlock(2, 2, alOp);
00173    }
00174    // 3D problem.
00175    else if(dim == 3)
00176    {
00177       lpA33 = getBlock(2, 2, alOp);
00178       lpB1  = getBlock(3, 0, alOp);
00179       lpB2  = getBlock(3, 1, alOp);
00180       lpB3  = getBlock(3, 2, alOp);
00181       lpB1t = getBlock(0, 3, alOp);
00182       lpB2t = getBlock(1, 3, alOp);
00183       lpB3t = getBlock(2, 3, alOp);
00184       lpC   = getBlock(3, 3, alOp);
00185    }
00186 
00187    // For problems using stabilized finite elements,
00188    // lpB1t, lpB2t and lpB3t are added linear operators. Extract original operators.
00189    LinearOp B1t = (rcp_dynamic_cast<const Thyra::DefaultAddedLinearOp<double> >(lpB1t))->getOp(0);
00190    LinearOp B2t = (rcp_dynamic_cast<const Thyra::DefaultAddedLinearOp<double> >(lpB2t))->getOp(0);
00191    LinearOp B3t;
00192    if(dim == 3)
00193    {
00194       B3t = (rcp_dynamic_cast<const Thyra::DefaultAddedLinearOp<double> >(lpB3t))->getOp(0);
00195    }
00196 
00197    //std::cout << Teuchos::describe(*lpC, Teuchos::VERB_EXTREME) << std::endl;
00198    // Check whether the finite elements are stable or not.
00199    state->isStabilized_ =(not isZeroOp(lpC));
00200    //state->isStabilized_ = false;
00201    //std::cout << state->isStabilized_ << std::endl;
00202 
00203    state->pressureMassMatrix_ = pressureMassMatrix_;
00204    // If pressure mass matrix is not set, use identity.
00205    if(state->pressureMassMatrix_ == Teuchos::null)
00206    {
00207       Teko_DEBUG_MSG("InvModALStrategy::initializeState(): Build identity type \""
00208             << getDiagonalName(scaleType_) << "\"", 1);
00209       state->invPressureMassMatrix_ = Thyra::identity<double>(lpB1->range());
00210    }
00211    // If the inverse of the pressure mass matrix is not set,
00212    // build it from the pressure mass matrix.
00213    else if(state->invPressureMassMatrix_ == Teuchos::null)
00214    {
00215       Teko_DEBUG_MSG("ModAL::initializeState(): Build Scaling <mass> type \""
00216             << getDiagonalName(scaleType_) << "\"", 1);
00217       state->invPressureMassMatrix_ = getInvDiagonalOp(pressureMassMatrix_, scaleType_);
00218    }
00219    // Else "invPressureMassMatrix_" should be set and there is no reason to rebuild it
00220    state->gamma_ = gamma_;
00221    //S_ = scale(1.0/gamma_, pressureMassMatrix_);
00222    std::cout << Teuchos::describe(*(state->invPressureMassMatrix_), Teuchos::VERB_EXTREME) << std::endl;
00223 
00224    // Build state variables: B_1^T*W^{-1}*B_1, A11p, etc.
00225    // B_1^T*W^{-1}*B_1 may not change so save it in the state.
00226    if(state->B1tMpB1_ == Teuchos::null)
00227       state->B1tMpB1_ = explicitMultiply(B1t, state->invPressureMassMatrix_, lpB1, state->B1tMpB1_);
00228    // Get the(1,1) block of the non-augmented matrix.
00229    // Recall alOp is augmented. So lpA11 = A11 + gamma B_1^T W^{-1} B_1.
00230    // Cast lpA11 as an added linear operator and get the first item.
00231    LinearOp A11 = (rcp_dynamic_cast<const Thyra::DefaultAddedLinearOp<double> >(lpA11))->getOp(0);
00232    state->A11p_ = explicitAdd(A11, scale(state->gamma_, state->B1tMpB1_), state->A11p_);
00233    //std::cout << Teuchos::describe(*(state->B1tMpB1_), Teuchos::VERB_EXTREME) << std::endl;
00234    Teko_DEBUG_MSG("Computed A11p", 10);
00235 
00236    if(state->B2tMpB2_ == Teuchos::null)
00237       state->B2tMpB2_ = explicitMultiply(B2t, state->invPressureMassMatrix_, lpB2, state->B2tMpB2_);
00238    LinearOp A22 = (rcp_dynamic_cast<const Thyra::DefaultAddedLinearOp<double> >(lpA22))->getOp(0);
00239    state->A22p_ = explicitAdd(A22, scale(state->gamma_, state->B2tMpB2_), state->A22p_);
00240    Teko_DEBUG_MSG("Computed A22p", 10);
00241 
00242    if(dim == 3)
00243    {
00244       if(state->B3tMpB3_ == Teuchos::null)
00245          state->B3tMpB3_ = explicitMultiply(B3t, state->invPressureMassMatrix_, lpB3, state->B3tMpB3_);
00246       LinearOp A33 = (rcp_dynamic_cast<const Thyra::DefaultAddedLinearOp<double> >(lpA33))->getOp(0);
00247       state->A33p_ = explicitAdd(A33, scale(state->gamma_, state->B3tMpB3_), state->A33p_);
00248       Teko_DEBUG_MSG("Computed A33p", 10);
00249    }
00250 
00251    // Inverse the Schur complement.
00252    if(state->isStabilized_)
00253    {
00254       if(state->S_ == Teuchos::null)
00255       {
00256          state->S_ = explicitAdd(scale(-1.0, lpC), scale(1.0/state->gamma_, pressureMassMatrix_), state->S_);
00257       }
00258       Teko_DEBUG_MSG("Computed S", 10);
00259    }
00260 
00261    state->setInitialized(true);
00262 }
00263 
00264 // Compute inverses.
00265 void InvModALStrategy::computeInverses(const BlockedLinearOp & alOp,
00266       ModALPrecondState *state) const
00267 {
00268    int dim = blockRowCount(alOp) - 1;
00269    TEUCHOS_ASSERT(dim == 2 || dim == 3);
00270 
00271    Teko_DEBUG_SCOPE("InvModALStrategy::computeInverses", 10);
00272    Teko_DEBUG_EXPR(Teuchos::Time invTimer(""));
00273 
00274    //(re)build the inverse of A11
00275    Teko_DEBUG_MSG("ModAL::computeInverses(): Building inv(A11)", 1);
00276    Teko_DEBUG_EXPR(invTimer.start(true));
00277 
00278    InverseLinearOp invA11p = state->getInverse("invA11p");
00279    if(invA11p == Teuchos::null)
00280    {
00281       invA11p = buildInverse(*invFactoryA_, state->A11p_);
00282       state->addInverse("invA11p", invA11p);
00283    }
00284    else
00285    {
00286       rebuildInverse(*invFactoryA_, state->A11p_, invA11p);
00287    }
00288 
00289    Teko_DEBUG_EXPR(invTimer.stop());
00290    Teko_DEBUG_MSG("ModAL::computeInverses GetInvA11 = " << invTimer.totalElapsedTime(), 1);
00291 
00292    //(re)build the inverse of A22
00293    Teko_DEBUG_MSG("ModAL::computeInverses(): Building inv(A22)", 2);
00294    Teko_DEBUG_EXPR(invTimer.start(true));
00295 
00296    InverseLinearOp invA22p = state->getInverse("invA22p");
00297    if(invA22p == Teuchos::null)
00298    {
00299       invA22p = buildInverse(*invFactoryA_, state->A22p_);
00300       state->addInverse("invA22p", invA22p);
00301    }
00302    else
00303    {
00304       rebuildInverse(*invFactoryA_, state->A22p_, invA22p);
00305    }
00306 
00307    Teko_DEBUG_EXPR(invTimer.stop());
00308    Teko_DEBUG_MSG("ModAL::computeInverses(): GetInvA22 = " << invTimer.totalElapsedTime(), 2);
00309 
00310    //(re)build the inverse of A33
00311    if(dim == 3)
00312    {
00313       Teko_DEBUG_MSG("ModAL::computeInverses Building inv(A33)", 3);
00314       Teko_DEBUG_EXPR(invTimer.start(true));
00315 
00316       InverseLinearOp invA33p = state->getInverse("invA33p");
00317       if(invA33p == Teuchos::null)
00318       {
00319          invA33p = buildInverse(*invFactoryA_, state->A33p_);
00320          state->addInverse("invA33p", invA33p);
00321       }
00322       else
00323       {
00324          rebuildInverse(*invFactoryA_, state->A33p_, invA33p);
00325       }
00326 
00327       Teko_DEBUG_EXPR(invTimer.stop());
00328       Teko_DEBUG_MSG("ModAL::computeInverses GetInvA33 = " << invTimer.totalElapsedTime(), 3);
00329    }
00330 
00331    //(re)build the inverse of S
00332    Teko_DEBUG_MSG("ModAL::computeInverses Building inv(S)", 4);
00333    Teko_DEBUG_EXPR(invTimer.start(true));
00334 
00335    // There are two ways to "invert" S.
00336    // The following method construct invS by InverseFactory invFactoryS_.
00337    // The other way is to use diagonal approximation,
00338    // which is done in ModALPreconditionerFactory.cpp.
00339    if(state->isStabilized_)
00340    {
00341       InverseLinearOp invS = state->getInverse("invS");
00342       if(invS == Teuchos::null)
00343       {
00344          invS = buildInverse(*invFactoryS_, state->S_);
00345          state->addInverse("invS", invS);
00346       }
00347       else
00348       {
00349          rebuildInverse(*invFactoryS_, state->S_, invS);
00350       }
00351    }
00352 
00353    Teko_DEBUG_EXPR(invTimer.stop());
00354    Teko_DEBUG_MSG("ModAL::computeInverses GetInvS = " << invTimer.totalElapsedTime(), 4);
00355 
00356 }
00357 
00358 } // end namespace NS
00359 
00360 } // end namespace Teko
 All Classes Files Functions Variables