7 #include "Thyra_DefaultDiagonalLinearOp.hpp" 8 #include "Thyra_EpetraThyraWrappers.hpp" 9 #include "Thyra_get_Epetra_Operator.hpp" 10 #include "Thyra_EpetraLinearOp.hpp" 11 #include "Thyra_DefaultAddedLinearOp_def.hpp" 12 #include "Thyra_DefaultIdentityLinearOp_decl.hpp" 14 #include "Epetra_Vector.h" 15 #include "Epetra_Map.h" 17 #include "EpetraExt_RowMatrixOut.h" 18 #include "EpetraExt_MultiVectorOut.h" 20 #include "Teuchos_Time.hpp" 24 #include "Teko_InvModALStrategy.hpp" 25 #include "Teko_ModALPreconditionerFactory.hpp" 28 using Teuchos::rcp_dynamic_cast;
29 using Teuchos::rcp_const_cast;
38 InvModALStrategy::InvModALStrategy() :
40 pressureMassMatrix_(
Teuchos::null), gamma_(0.05),
41 scaleType_(
Diagonal), isSymmetric_(true)
46 InvModALStrategy::InvModALStrategy(
const Teuchos::RCP<InverseFactory> & factory) :
47 invFactoryA_(factory), invFactoryS_(factory),
48 pressureMassMatrix_(
Teuchos::null), gamma_(0.05),
49 scaleType_(
Diagonal), isSymmetric_(true)
54 InvModALStrategy::InvModALStrategy(
const Teuchos::RCP<InverseFactory> & invFactoryA,
55 const Teuchos::RCP<InverseFactory> & invFactoryS) :
56 invFactoryA_(invFactoryA), invFactoryS_(invFactoryS),
57 pressureMassMatrix_(
Teuchos::null), gamma_(0.05),
58 scaleType_(
Diagonal), isSymmetric_(true)
63 InvModALStrategy::InvModALStrategy(
const Teuchos::RCP<InverseFactory> & invFactory,
64 LinearOp & pressureMassMatrix) :
65 invFactoryA_(invFactory), invFactoryS_(invFactory),
66 pressureMassMatrix_(pressureMassMatrix), gamma_(0.05),
67 scaleType_(
Diagonal), isSymmetric_(true)
72 InvModALStrategy::InvModALStrategy(
const Teuchos::RCP<InverseFactory> & invFactoryA,
73 const Teuchos::RCP<InverseFactory> & invFactoryS, LinearOp & pressureMassMatrix) :
74 invFactoryA_(invFactoryA), invFactoryS_(invFactoryS),
75 pressureMassMatrix_(pressureMassMatrix), gamma_(0.05),
76 scaleType_(
Diagonal), isSymmetric_(true)
81 LinearOp InvModALStrategy::getInvA11p(BlockPreconditionerState & state)
const 83 return state.getInverse(
"invA11p");
86 LinearOp InvModALStrategy::getInvA22p(BlockPreconditionerState & state)
const 88 return state.getInverse(
"invA22p");
91 LinearOp InvModALStrategy::getInvA33p(BlockPreconditionerState & state)
const 93 return state.getInverse(
"invA33p");
96 LinearOp InvModALStrategy::getInvS(BlockPreconditionerState & state)
const 98 return state.getInverse(
"invS");
102 void InvModALStrategy::setPressureMassMatrix(
const LinearOp & pressureMassMatrix)
104 pressureMassMatrix_ = pressureMassMatrix;
108 void InvModALStrategy::setGamma(
double gamma)
110 TEUCHOS_ASSERT(gamma > 0.0);
114 void InvModALStrategy::buildState(
const BlockedLinearOp & alOp,
115 BlockPreconditionerState & state)
const 117 Teko_DEBUG_SCOPE(
"InvModALStrategy::buildState", 10);
119 ModALPrecondState * modALState =
dynamic_cast<ModALPrecondState*
>(&state);
120 TEUCHOS_ASSERT(modALState != NULL);
123 if(not modALState->isInitialized())
125 Teko_DEBUG_EXPR(Teuchos::Time timer(
""));
129 Teko_DEBUG_SCOPE(
"ModAL::buildState: Initializing state object", 1);
130 Teko_DEBUG_EXPR(timer.start(
true));
134 Teko_DEBUG_EXPR(timer.stop());
135 Teko_DEBUG_MSG(
"ModAL::buildState: BuildOpsTime = " << timer.totalElapsedTime(), 1);
140 Teko_DEBUG_SCOPE(
"ModAL::buildState: Computing inverses", 1);
141 Teko_DEBUG_EXPR(timer.start(
true));
145 Teko_DEBUG_EXPR(timer.stop());
146 Teko_DEBUG_MSG(
"ModAL::buildState: BuildInvTime = " << timer.totalElapsedTime(), 1);
152 void InvModALStrategy::initializeState(
const BlockedLinearOp & alOp,
153 ModALPrecondState *state)
const 155 Teko_DEBUG_SCOPE(
"InvModALStrategy::initializeState", 10);
159 TEUCHOS_ASSERT(dim == 2 || dim == 3);
161 LinearOp lpA11 =
getBlock(0, 0, alOp);
162 LinearOp lpA22 =
getBlock(1, 1, alOp);
163 LinearOp lpA33, lpB1, lpB2, lpB3, lpB1t, lpB2t, lpB3t, lpC;
189 LinearOp B1t = (rcp_dynamic_cast<
const Thyra::DefaultAddedLinearOp<double> >(lpB1t))->getOp(0);
190 LinearOp B2t = (rcp_dynamic_cast<
const Thyra::DefaultAddedLinearOp<double> >(lpB2t))->getOp(0);
194 B3t = (rcp_dynamic_cast<
const Thyra::DefaultAddedLinearOp<double> >(lpB3t))->getOp(0);
199 state->isStabilized_ =(not isZeroOp(lpC));
203 state->pressureMassMatrix_ = pressureMassMatrix_;
205 if(state->pressureMassMatrix_ == Teuchos::null)
207 Teko_DEBUG_MSG(
"InvModALStrategy::initializeState(): Build identity type \"" 208 << getDiagonalName(scaleType_) <<
"\"", 1);
209 state->invPressureMassMatrix_ = Thyra::identity<double>(lpB1->range());
213 else if(state->invPressureMassMatrix_ == Teuchos::null)
215 Teko_DEBUG_MSG(
"ModAL::initializeState(): Build Scaling <mass> type \"" 216 << getDiagonalName(scaleType_) <<
"\"", 1);
217 state->invPressureMassMatrix_ = getInvDiagonalOp(pressureMassMatrix_, scaleType_);
220 state->gamma_ = gamma_;
222 std::cout << Teuchos::describe(*(state->invPressureMassMatrix_), Teuchos::VERB_EXTREME) << std::endl;
226 if(state->B1tMpB1_ == Teuchos::null)
227 state->B1tMpB1_ = explicitMultiply(B1t, state->invPressureMassMatrix_, lpB1, state->B1tMpB1_);
231 LinearOp A11 = (rcp_dynamic_cast<
const Thyra::DefaultAddedLinearOp<double> >(lpA11))->getOp(0);
232 state->A11p_ = explicitAdd(A11,
scale(state->gamma_, state->B1tMpB1_), state->A11p_);
234 Teko_DEBUG_MSG(
"Computed A11p", 10);
236 if(state->B2tMpB2_ == Teuchos::null)
237 state->B2tMpB2_ = explicitMultiply(B2t, state->invPressureMassMatrix_, lpB2, state->B2tMpB2_);
238 LinearOp A22 = (rcp_dynamic_cast<
const Thyra::DefaultAddedLinearOp<double> >(lpA22))->getOp(0);
239 state->A22p_ = explicitAdd(A22,
scale(state->gamma_, state->B2tMpB2_), state->A22p_);
240 Teko_DEBUG_MSG(
"Computed A22p", 10);
244 if(state->B3tMpB3_ == Teuchos::null)
245 state->B3tMpB3_ = explicitMultiply(B3t, state->invPressureMassMatrix_, lpB3, state->B3tMpB3_);
246 LinearOp A33 = (rcp_dynamic_cast<
const Thyra::DefaultAddedLinearOp<double> >(lpA33))->getOp(0);
247 state->A33p_ = explicitAdd(A33,
scale(state->gamma_, state->B3tMpB3_), state->A33p_);
248 Teko_DEBUG_MSG(
"Computed A33p", 10);
252 if(state->isStabilized_)
254 if(state->S_ == Teuchos::null)
256 state->S_ = explicitAdd(
scale(-1.0, lpC),
scale(1.0/state->gamma_, pressureMassMatrix_), state->S_);
258 Teko_DEBUG_MSG(
"Computed S", 10);
261 state->setInitialized(
true);
265 void InvModALStrategy::computeInverses(
const BlockedLinearOp & alOp,
266 ModALPrecondState *state)
const 269 TEUCHOS_ASSERT(dim == 2 || dim == 3);
271 Teko_DEBUG_SCOPE(
"InvModALStrategy::computeInverses", 10);
272 Teko_DEBUG_EXPR(Teuchos::Time invTimer(
""));
275 Teko_DEBUG_MSG(
"ModAL::computeInverses(): Building inv(A11)", 1);
276 Teko_DEBUG_EXPR(invTimer.start(
true));
278 InverseLinearOp invA11p = state->getInverse(
"invA11p");
279 if(invA11p == Teuchos::null)
282 state->addInverse(
"invA11p", invA11p);
289 Teko_DEBUG_EXPR(invTimer.stop());
290 Teko_DEBUG_MSG(
"ModAL::computeInverses GetInvA11 = " << invTimer.totalElapsedTime(), 1);
293 Teko_DEBUG_MSG(
"ModAL::computeInverses(): Building inv(A22)", 2);
294 Teko_DEBUG_EXPR(invTimer.start(
true));
296 InverseLinearOp invA22p = state->getInverse(
"invA22p");
297 if(invA22p == Teuchos::null)
300 state->addInverse(
"invA22p", invA22p);
307 Teko_DEBUG_EXPR(invTimer.stop());
308 Teko_DEBUG_MSG(
"ModAL::computeInverses(): GetInvA22 = " << invTimer.totalElapsedTime(), 2);
313 Teko_DEBUG_MSG(
"ModAL::computeInverses Building inv(A33)", 3);
314 Teko_DEBUG_EXPR(invTimer.start(
true));
316 InverseLinearOp invA33p = state->getInverse(
"invA33p");
317 if(invA33p == Teuchos::null)
320 state->addInverse(
"invA33p", invA33p);
327 Teko_DEBUG_EXPR(invTimer.stop());
328 Teko_DEBUG_MSG(
"ModAL::computeInverses GetInvA33 = " << invTimer.totalElapsedTime(), 3);
332 Teko_DEBUG_MSG(
"ModAL::computeInverses Building inv(S)", 4);
333 Teko_DEBUG_EXPR(invTimer.start(
true));
339 if(state->isStabilized_)
341 InverseLinearOp invS = state->getInverse(
"invS");
342 if(invS == Teuchos::null)
345 state->addInverse(
"invS", invS);
353 Teko_DEBUG_EXPR(invTimer.stop());
354 Teko_DEBUG_MSG(
"ModAL::computeInverses GetInvS = " << invTimer.totalElapsedTime(), 4);
Specifies that just the diagonal is used.
int blockRowCount(const BlockedLinearOp &blo)
Get the row count in a block linear operator.
InverseLinearOp buildInverse(const InverseFactory &factory, const LinearOp &A)
Build an inverse operator using a factory and a linear operator.
void scale(const double alpha, MultiVector &x)
Scale a multivector by a constant.
void rebuildInverse(const InverseFactory &factory, const LinearOp &A, InverseLinearOp &invA)
virtual void initializeState(const BlockedLinearOp &A, LSCPrecondState *state) const
Initialize the state object using this blocked linear operator.
void computeInverses(const BlockedLinearOp &A, LSCPrecondState *state) const
MultiVector getBlock(int i, const BlockedMultiVector &bmv)
Get the ith block from a BlockedMultiVector object.