Teko Version of the Day
Loading...
Searching...
No Matches
Teko_ALOperator.cpp
1/*
2 * Author: Zhen Wang
3 * Email: wangz@ornl.gov
4 * zhen.wang@alum.emory.edu
5 */
6
7#include "Teko_BlockedEpetraOperator.hpp"
8#include "Teko_BlockedMappingStrategy.hpp"
9#include "Teko_ReorderedMappingStrategy.hpp"
10
11#include "Teuchos_VerboseObject.hpp"
12
13#include "Thyra_LinearOpBase.hpp"
14#include "Thyra_EpetraLinearOp.hpp"
15#include "Thyra_EpetraThyraWrappers.hpp"
16#include "Thyra_DefaultProductMultiVector.hpp"
17#include "Thyra_DefaultProductVectorSpace.hpp"
18#include "Thyra_DefaultBlockedLinearOp.hpp"
19#include "Thyra_get_Epetra_Operator.hpp"
20
21#include "EpetraExt_MultiVectorOut.h"
22#include "EpetraExt_RowMatrixOut.h"
23
24#include "Teko_Utilities.hpp"
25
26#include "Teko_ALOperator.hpp"
27
28namespace Teko
29{
30
31namespace NS
32{
33
34ALOperator::ALOperator(const std::vector<std::vector<int> > & vars,
35 const Teuchos::RCP<Epetra_Operator> & content,
36 LinearOp pressureMassMatrix, double gamma,
37 const std::string & label) :
38 Teko::Epetra::BlockedEpetraOperator(vars, content, label),
39 pressureMassMatrix_(pressureMassMatrix), gamma_(gamma)
40{
41 checkDim(vars);
42 SetContent(vars, content);
44}
45
46ALOperator::ALOperator(const std::vector<std::vector<int> > & vars,
47 const Teuchos::RCP<Epetra_Operator> & content,
48 double gamma, const std::string & label) :
49 Teko::Epetra::BlockedEpetraOperator(vars, content, label),
50 pressureMassMatrix_(Teuchos::null), gamma_(gamma)
51{
52 checkDim(vars);
53 SetContent(vars, content);
55}
56
57void
58ALOperator::setPressureMassMatrix(LinearOp pressureMassMatrix)
59{
60 if(pressureMassMatrix != Teuchos::null)
61 {
62 pressureMassMatrix_ = pressureMassMatrix;
63 }
64}
65
66void
68{
69 TEUCHOS_ASSERT(gamma > 0.0);
70 gamma_ = gamma;
71}
72
73const Teuchos::RCP<const Epetra_Operator>
74ALOperator::GetBlock(int i, int j) const
75{
76 const Teuchos::RCP<Thyra::BlockedLinearOpBase<double> > blkOp
77 = Teuchos::rcp_dynamic_cast< Thyra::BlockedLinearOpBase<double> >
78 (blockedOperator_, true);
79
80 return Thyra::get_Epetra_Operator(*blkOp->getBlock(i, j));
81}
82
83void
84ALOperator::checkDim(const std::vector<std::vector<int> > & vars)
85{
86 dim_ = vars.size() - 1;
87 TEUCHOS_ASSERT(dim_ == 2 || dim_ == 3);
88}
89
90void
92{
93 TEUCHOS_ASSERT(blockedMapping_ != Teuchos::null);
94
95 // Get an Epetra_CrsMatrix.
96 const Teuchos::RCP<const Epetra_CrsMatrix> crsContent
97 = Teuchos::rcp_dynamic_cast<const Epetra_CrsMatrix>(fullContent_);
98
99 // Ask the strategy to build the Thyra operator for you.
100 if(blockedOperator_ == Teuchos::null)
101 {
102 blockedOperator_
103 = blockedMapping_->buildBlockedThyraOp(crsContent, label_);
104 }
105 else
106 {
107 const Teuchos::RCP<Thyra::BlockedLinearOpBase<double> > blkOp
108 = Teuchos::rcp_dynamic_cast<Thyra::BlockedLinearOpBase<double> >
109 (blockedOperator_, true);
110 blockedMapping_->rebuildBlockedThyraOp(crsContent, blkOp);
111 }
112
113 // Extract blocks.
114 const Teuchos::RCP<Thyra::BlockedLinearOpBase<double> > blkOp
115 = Teuchos::rcp_dynamic_cast<Thyra::BlockedLinearOpBase<double> >
116 (blockedOperator_, true);
117 numBlockRows_ = blkOp->productRange()->numBlocks();
118 Teuchos::RCP<const Thyra::LinearOpBase<double> > blockedOpBlocks[4][4];
119 for(int i = 0; i <= dim_; i++)
120 {
121 for(int j = 0; j <= dim_; j++)
122 {
123 blockedOpBlocks[i][j] = blkOp->getBlock(i, j);
124 }
125 }
126
127 // Pressure mass matrix.
128 if(pressureMassMatrix_ != Teuchos::null)
129 {
131 }
132 // We need the size of the sub-block to build the identity matrix.
133 else
134 {
135 std::cout << "Pressure mass matrix is null. Use identity." << std::endl;
137 = Thyra::identity<double>(blockedOpBlocks[dim_][0]->range());
139 = Thyra::identity<double>(blockedOpBlocks[dim_][0]->range());
140 }
141
142 // Build the AL operator.
143 Teuchos::RCP<Thyra::DefaultBlockedLinearOp<double> > alOperator
144 = Thyra::defaultBlockedLinearOp<double>();
145 alOperator->beginBlockFill(dim_ + 1, dim_ + 1);
146
147 // Set blocks for the velocity parts and gradient.
148 for(int i = 0; i < dim_; i++)
149 {
150 for(int j = 0; j < dim_; j++)
151 {
152 // build the blocks and place it the right location
153 alOperator->setBlock(i, j,
154 Thyra::add(blockedOpBlocks[i][j],
155 Thyra::scale(gamma_,
156 Thyra::multiply(blockedOpBlocks[i][dim_],
157 invPressureMassMatrix_, blockedOpBlocks[dim_][j]))));
158 } // end for j
159 } // end for i
160
161 // Last row. Divergence and (possible) stabilization matrix.
162 for(int j = 0; j <= dim_; j++)
163 {
164 alOperator->setBlock(dim_, j, blockedOpBlocks[dim_][j]);
165 }
166
167 // Last column. Gradient.
168 for(int i = 0; i < dim_; i++)
169 {
170 alOperator->setBlock(i, dim_,
171 Thyra::add(blockedOpBlocks[i][dim_],
172 Thyra::scale(gamma_,
173 Thyra::multiply(blockedOpBlocks[i][dim_],
174 invPressureMassMatrix_,blockedOpBlocks[dim_][dim_]))));
175 }
176
177 alOperator->endBlockFill();
178
179 // Set whatever is returned.
180 SetOperator(alOperator, false);
181
182 // Set operator for augmenting the right-hand side.
183 Teuchos::RCP<Thyra::DefaultBlockedLinearOp<double> > alOpRhs_
184 = Thyra::defaultBlockedLinearOp<double>();
185 alOpRhs_->beginBlockFill(dim_ + 1, dim_ + 1);
186
187 // Identity matrices on the main diagonal.
188 for(int i = 0; i < dim_; i++)
189 {
190 // build the blocks and place it the right location
191 alOpRhs_->setBlock(i, i,
192 Thyra::identity<double>(blockedOpBlocks[0][0]->range()));
193 } // end for i
194 alOpRhs_->setBlock(dim_, dim_,
195 Thyra::identity<double>(blockedOpBlocks[dim_][dim_]->range()));
196
197 // Last column.
198 for(int i = 0; i < dim_; i++)
199 {
200 alOpRhs_->setBlock(i, dim_,
201 Thyra::scale(gamma_,
202 Thyra::multiply(blockedOpBlocks[i][dim_], invPressureMassMatrix_)));
203 }
204
205 alOpRhs_->endBlockFill();
206
207 alOperatorRhs_ = alOpRhs_;
208
209 // reorder if necessary
210 if(reorderManager_ != Teuchos::null)
211 Reorder(*reorderManager_);
212}
213
214void
215ALOperator::augmentRHS(const Epetra_MultiVector & b, Epetra_MultiVector & bAugmented)
216{
217 Teuchos::RCP<const Teko::Epetra::MappingStrategy> mapping
218 = this->getMapStrategy();
219 Teuchos::RCP<Thyra::MultiVectorBase<double> > bThyra
220 = Thyra::createMembers(thyraOp_->range(), b.NumVectors());
221 Teuchos::RCP<Thyra::MultiVectorBase<double> > bThyraAugmented
222 = Thyra::createMembers(thyraOp_->range(), b.NumVectors());
223 //std::cout << Teuchos::describe(*bThyra, Teuchos::VERB_EXTREME) << std::endl;
224 // Copy Epetra vector to Thyra vector.
225 mapping->copyEpetraIntoThyra(b, bThyra.ptr());
226 // Apply operator.
227 alOperatorRhs_->apply(Thyra::NOTRANS, *bThyra, bThyraAugmented.ptr(), 1.0, 0.0);
228 // Copy Thyra vector to Epetra vector.
229 mapping->copyThyraIntoEpetra(bThyraAugmented, bAugmented);
230}
231
232} // end namespace NS
233
234} // end namespace Teko
virtual void SetContent(const std::vector< std::vector< int > > &vars, const Teuchos::RCP< const Epetra_Operator > &content)
void Reorder(const BlockReorderManager &brm)
const RCP< const MappingStrategy > getMapStrategy() const
Get the mapping strategy for this wrapper (translate between Thyra and Epetra)
void setPressureMassMatrix(LinearOp pressureMassMatrix)
Teuchos::RCP< Thyra::LinearOpBase< double > > alOperatorRhs_
const Teuchos::RCP< const Epetra_Operator > GetBlock(int i, int j) const
void setGamma(double gamma)
void checkDim(const std::vector< std::vector< int > > &vars)
void augmentRHS(const Epetra_MultiVector &b, Epetra_MultiVector &bAugmented)
ALOperator(const std::vector< std::vector< int > > &vars, const Teuchos::RCP< Epetra_Operator > &content, LinearOp pressureMassMatrix, double gamma=0.05, const std::string &label="<ANYM>")