Teko Version of the Day
Teko_ALOperator.cpp
00001 /*
00002  * Author: Zhen Wang
00003  * Email: wangz@ornl.gov
00004  *        zhen.wang@alum.emory.edu
00005  */
00006 
00007 #include "Teko_BlockedEpetraOperator.hpp"
00008 #include "Teko_BlockedMappingStrategy.hpp"
00009 #include "Teko_ReorderedMappingStrategy.hpp"
00010 
00011 #include "Teuchos_VerboseObject.hpp"
00012 
00013 #include "Thyra_LinearOpBase.hpp"
00014 #include "Thyra_EpetraLinearOp.hpp"
00015 #include "Thyra_EpetraThyraWrappers.hpp"
00016 #include "Thyra_DefaultProductMultiVector.hpp"
00017 #include "Thyra_DefaultProductVectorSpace.hpp"
00018 #include "Thyra_DefaultBlockedLinearOp.hpp"
00019 #include "Thyra_get_Epetra_Operator.hpp"
00020 
00021 #include "EpetraExt_MultiVectorOut.h"
00022 #include "EpetraExt_RowMatrixOut.h"
00023 
00024 #include "Teko_Utilities.hpp"
00025 
00026 #include "Teko_ALOperator.hpp"
00027 
00028 namespace Teko
00029 {
00030 
00031 namespace NS
00032 {
00033 
00034 ALOperator::ALOperator(const std::vector<std::vector<int> > & vars,
00035       const Teuchos::RCP<Epetra_Operator> & content,
00036       LinearOp pressureMassMatrix, double gamma,
00037       const std::string & label) :
00038       Teko::Epetra::BlockedEpetraOperator(vars, content, label),
00039       pressureMassMatrix_(pressureMassMatrix), gamma_(gamma)
00040 {
00041    checkDim(vars);
00042    SetContent(vars, content);
00043    BuildALOperator();
00044 }
00045 
00046 ALOperator::ALOperator(const std::vector<std::vector<int> > & vars,
00047       const Teuchos::RCP<Epetra_Operator> & content,
00048       double gamma, const std::string & label) :
00049       Teko::Epetra::BlockedEpetraOperator(vars, content, label),
00050       pressureMassMatrix_(Teuchos::null), gamma_(gamma)
00051 {
00052    checkDim(vars);
00053    SetContent(vars, content);
00054    BuildALOperator();
00055 }
00056 
00057 void
00058 ALOperator::setPressureMassMatrix(LinearOp pressureMassMatrix)
00059 {
00060    if(pressureMassMatrix != Teuchos::null)
00061    {
00062       pressureMassMatrix_ = pressureMassMatrix;
00063    }
00064 }
00065 
00066 void
00067 ALOperator::setGamma(double gamma)
00068 {
00069    TEUCHOS_ASSERT(gamma > 0.0);
00070    gamma_ = gamma;
00071 }
00072 
00073 const Teuchos::RCP<const Epetra_Operator>
00074 ALOperator::GetBlock(int i, int j) const
00075 {
00076    const Teuchos::RCP<Thyra::BlockedLinearOpBase<double> > blkOp
00077          = Teuchos::rcp_dynamic_cast< Thyra::BlockedLinearOpBase<double> >
00078          (blockedOperator_, true);
00079 
00080    return Thyra::get_Epetra_Operator(*blkOp->getBlock(i, j));
00081 }
00082 
00083 void
00084 ALOperator::checkDim(const std::vector<std::vector<int> > & vars)
00085 {
00086    dim_ = vars.size() - 1;
00087    TEUCHOS_ASSERT(dim_ == 2 || dim_ == 3);
00088 }
00089 
00090 void
00091 ALOperator::BuildALOperator()
00092 {
00093    TEUCHOS_ASSERT(blockedMapping_ != Teuchos::null);
00094 
00095    // Get an Epetra_CrsMatrix.
00096    const Teuchos::RCP<const Epetra_CrsMatrix> crsContent
00097       = Teuchos::rcp_dynamic_cast<const Epetra_CrsMatrix>(fullContent_);
00098 
00099    // Ask the strategy to build the Thyra operator for you.
00100    if(blockedOperator_ == Teuchos::null)
00101    {
00102       blockedOperator_
00103          = blockedMapping_->buildBlockedThyraOp(crsContent, label_);
00104    }
00105    else
00106    {
00107       const Teuchos::RCP<Thyra::BlockedLinearOpBase<double> > blkOp
00108             = Teuchos::rcp_dynamic_cast<Thyra::BlockedLinearOpBase<double> >
00109             (blockedOperator_, true);
00110       blockedMapping_->rebuildBlockedThyraOp(crsContent, blkOp);
00111    }
00112 
00113    // Extract blocks.
00114    const Teuchos::RCP<Thyra::BlockedLinearOpBase<double> > blkOp
00115       = Teuchos::rcp_dynamic_cast<Thyra::BlockedLinearOpBase<double> >
00116       (blockedOperator_, true);
00117    numBlockRows_ = blkOp->productRange()->numBlocks();
00118    Teuchos::RCP<const Thyra::LinearOpBase<double> > blockedOpBlocks[4][4];
00119    for(int i = 0; i <= dim_; i++)
00120    {
00121       for(int j = 0; j <= dim_; j++)
00122       {
00123          blockedOpBlocks[i][j] = blkOp->getBlock(i, j);
00124       }
00125    }
00126 
00127    // Pressure mass matrix.
00128    if(pressureMassMatrix_ != Teuchos::null)
00129    {
00130       invPressureMassMatrix_ = getInvDiagonalOp(pressureMassMatrix_);
00131    }
00132    // We need the size of the sub-block to build the identity matrix.
00133    else
00134    {
00135       std::cout << "Pressure mass matrix is null. Use identity." << std::endl;
00136       pressureMassMatrix_
00137             = Thyra::identity<double>(blockedOpBlocks[dim_][0]->range());
00138       invPressureMassMatrix_
00139             = Thyra::identity<double>(blockedOpBlocks[dim_][0]->range());
00140    }
00141 
00142    // Build the AL operator.
00143    Teuchos::RCP<Thyra::DefaultBlockedLinearOp<double> > alOperator_
00144       = Thyra::defaultBlockedLinearOp<double>();
00145    alOperator_->beginBlockFill(dim_ + 1, dim_ + 1);
00146 
00147    // Set blocks for the velocity parts and gradient.
00148    for(int i = 0; i < dim_; i++)
00149    {
00150       for(int j = 0; j < dim_; j++)
00151       {
00152          // build the blocks and place it the right location
00153          alOperator_->setBlock(i, j,
00154                Thyra::add(blockedOpBlocks[i][j],
00155                Thyra::scale(gamma_,
00156                Thyra::multiply(blockedOpBlocks[i][dim_],
00157                invPressureMassMatrix_, blockedOpBlocks[dim_][j]))));
00158       } // end for j
00159    } // end for i
00160 
00161    // Last row. Divergence and (possible) stabilization matrix.
00162    for(int j = 0; j <= dim_; j++)
00163    {
00164       alOperator_->setBlock(dim_, j, blockedOpBlocks[dim_][j]);
00165    }
00166 
00167    // Last column. Gradient.
00168    for(int i = 0; i < dim_; i++)
00169    {
00170       alOperator_->setBlock(i, dim_,
00171             Thyra::add(blockedOpBlocks[i][dim_],
00172             Thyra::scale(gamma_,
00173             Thyra::multiply(blockedOpBlocks[i][dim_],
00174             invPressureMassMatrix_,blockedOpBlocks[dim_][dim_]))));
00175    }
00176 
00177    alOperator_->endBlockFill();
00178 
00179    // Set whatever is returned.
00180    SetOperator(alOperator_, false);
00181 
00182    // Set operator for augmenting the right-hand side.
00183    Teuchos::RCP<Thyra::DefaultBlockedLinearOp<double> > alOpRhs_
00184          = Thyra::defaultBlockedLinearOp<double>();
00185    alOpRhs_->beginBlockFill(dim_ + 1, dim_ + 1);
00186 
00187    // Identity matrices on the main diagonal.
00188    for(int i = 0; i < dim_; i++)
00189    {
00190       // build the blocks and place it the right location
00191       alOpRhs_->setBlock(i, i,
00192             Thyra::identity<double>(blockedOpBlocks[0][0]->range()));
00193    } // end for i
00194    alOpRhs_->setBlock(dim_, dim_,
00195          Thyra::identity<double>(blockedOpBlocks[dim_][dim_]->range()));
00196 
00197    // Last column.
00198    for(int i = 0; i < dim_; i++)
00199    {
00200       alOpRhs_->setBlock(i, dim_,
00201             Thyra::scale(gamma_,
00202             Thyra::multiply(blockedOpBlocks[i][dim_], invPressureMassMatrix_)));
00203    }
00204 
00205    alOpRhs_->endBlockFill();
00206 
00207    alOperatorRhs_ = alOpRhs_;
00208 
00209    // reorder if necessary
00210    if(reorderManager_ != Teuchos::null)
00211       Reorder(*reorderManager_);
00212 }
00213 
00214 void
00215 ALOperator::augmentRHS(const Epetra_MultiVector & b, Epetra_MultiVector & bAugmented)
00216 {
00217    Teuchos::RCP<const Teko::Epetra::MappingStrategy> mapping
00218          = this->getMapStrategy();
00219    Teuchos::RCP<Thyra::MultiVectorBase<double> > bThyra
00220          = Thyra::createMembers(thyraOp_->range(), b.NumVectors());
00221    Teuchos::RCP<Thyra::MultiVectorBase<double> > bThyraAugmented
00222          = Thyra::createMembers(thyraOp_->range(), b.NumVectors());
00223    //std::cout << Teuchos::describe(*bThyra, Teuchos::VERB_EXTREME) << std::endl;
00224    // Copy Epetra vector to Thyra vector.
00225    mapping->copyEpetraIntoThyra(b, bThyra.ptr());
00226    // Apply operator.
00227    alOperatorRhs_->apply(Thyra::NOTRANS, *bThyra, bThyraAugmented.ptr(), 1.0, 0.0);
00228    // Copy Thyra vector to Epetra vector.
00229    mapping->copyThyraIntoEpetra(bThyraAugmented, bAugmented);
00230 }
00231 
00232 } // end namespace NS
00233 
00234 } // end namespace Teko
 All Classes Files Functions Variables