7#include "Thyra_DefaultDiagonalLinearOp.hpp"
8#include "Thyra_DefaultAddedLinearOp_def.hpp"
9#include "Thyra_DefaultIdentityLinearOp_decl.hpp"
11#include "Teuchos_Time.hpp"
15#include "Teko_InvModALStrategy.hpp"
16#include "Teko_ModALPreconditionerFactory.hpp"
19using Teuchos::rcp_dynamic_cast;
20using Teuchos::rcp_const_cast;
29InvModALStrategy::InvModALStrategy() :
30 invFactoryA_(Teuchos::null), invFactoryS_(Teuchos::null),
31 pressureMassMatrix_(Teuchos::null), gamma_(0.05),
32 scaleType_(
Diagonal), isSymmetric_(true)
37InvModALStrategy::InvModALStrategy(
const Teuchos::RCP<InverseFactory> & factory) :
38 invFactoryA_(factory), invFactoryS_(factory),
39 pressureMassMatrix_(Teuchos::null), gamma_(0.05),
40 scaleType_(
Diagonal), isSymmetric_(true)
45InvModALStrategy::InvModALStrategy(
const Teuchos::RCP<InverseFactory> & invFactoryA,
46 const Teuchos::RCP<InverseFactory> & invFactoryS) :
47 invFactoryA_(invFactoryA), invFactoryS_(invFactoryS),
48 pressureMassMatrix_(Teuchos::null), gamma_(0.05),
49 scaleType_(
Diagonal), isSymmetric_(true)
54InvModALStrategy::InvModALStrategy(
const Teuchos::RCP<InverseFactory> & invFactory,
55 LinearOp & pressureMassMatrix) :
56 invFactoryA_(invFactory), invFactoryS_(invFactory),
57 pressureMassMatrix_(pressureMassMatrix), gamma_(0.05),
58 scaleType_(
Diagonal), isSymmetric_(true)
63InvModALStrategy::InvModALStrategy(
const Teuchos::RCP<InverseFactory> & invFactoryA,
64 const Teuchos::RCP<InverseFactory> & invFactoryS, LinearOp & pressureMassMatrix) :
65 invFactoryA_(invFactoryA), invFactoryS_(invFactoryS),
66 pressureMassMatrix_(pressureMassMatrix), gamma_(0.05),
67 scaleType_(
Diagonal), isSymmetric_(true)
72LinearOp InvModALStrategy::getInvA11p(BlockPreconditionerState & state)
const
74 return state.getInverse(
"invA11p");
77LinearOp InvModALStrategy::getInvA22p(BlockPreconditionerState & state)
const
79 return state.getInverse(
"invA22p");
82LinearOp InvModALStrategy::getInvA33p(BlockPreconditionerState & state)
const
84 return state.getInverse(
"invA33p");
87LinearOp InvModALStrategy::getInvS(BlockPreconditionerState & state)
const
89 return state.getInverse(
"invS");
93void InvModALStrategy::setPressureMassMatrix(
const LinearOp & pressureMassMatrix)
95 pressureMassMatrix_ = pressureMassMatrix;
99void InvModALStrategy::setGamma(
double gamma)
101 TEUCHOS_ASSERT(gamma > 0.0);
105void InvModALStrategy::buildState(
const BlockedLinearOp & alOp,
106 BlockPreconditionerState & state)
const
108 Teko_DEBUG_SCOPE(
"InvModALStrategy::buildState", 10);
110 ModALPrecondState * modALState =
dynamic_cast<ModALPrecondState*
>(&state);
111 TEUCHOS_ASSERT(modALState != NULL);
114 if(not modALState->isInitialized())
116 Teko_DEBUG_EXPR(Teuchos::Time timer(
""));
120 Teko_DEBUG_SCOPE(
"ModAL::buildState: Initializing state object", 1);
121 Teko_DEBUG_EXPR(timer.start(
true));
123 initializeState(alOp, modALState);
125 Teko_DEBUG_EXPR(timer.stop());
126 Teko_DEBUG_MSG(
"ModAL::buildState: BuildOpsTime = " << timer.totalElapsedTime(), 1);
131 Teko_DEBUG_SCOPE(
"ModAL::buildState: Computing inverses", 1);
132 Teko_DEBUG_EXPR(timer.start(
true));
134 computeInverses(alOp, modALState);
136 Teko_DEBUG_EXPR(timer.stop());
137 Teko_DEBUG_MSG(
"ModAL::buildState: BuildInvTime = " << timer.totalElapsedTime(), 1);
143void InvModALStrategy::initializeState(
const BlockedLinearOp & alOp,
144 ModALPrecondState *state)
const
146 Teko_DEBUG_SCOPE(
"InvModALStrategy::initializeState", 10);
150 TEUCHOS_ASSERT(dim == 2 || dim == 3);
152 LinearOp lpA11 =
getBlock(0, 0, alOp);
153 LinearOp lpA22 =
getBlock(1, 1, alOp);
154 LinearOp lpA33, lpB1, lpB2, lpB3, lpB1t, lpB2t, lpB3t, lpC;
180 LinearOp B1t = (rcp_dynamic_cast<const Thyra::DefaultAddedLinearOp<double> >(lpB1t))->getOp(0);
181 LinearOp B2t = (rcp_dynamic_cast<const Thyra::DefaultAddedLinearOp<double> >(lpB2t))->getOp(0);
185 B3t = (rcp_dynamic_cast<const Thyra::DefaultAddedLinearOp<double> >(lpB3t))->getOp(0);
190 state->isStabilized_ =(not isZeroOp(lpC));
194 state->pressureMassMatrix_ = pressureMassMatrix_;
196 if(state->pressureMassMatrix_ == Teuchos::null)
198 Teko_DEBUG_MSG(
"InvModALStrategy::initializeState(): Build identity type \""
199 << getDiagonalName(scaleType_) <<
"\"", 1);
200 state->invPressureMassMatrix_ = Thyra::identity<double>(lpB1->range());
204 else if(state->invPressureMassMatrix_ == Teuchos::null)
206 Teko_DEBUG_MSG(
"ModAL::initializeState(): Build Scaling <mass> type \""
207 << getDiagonalName(scaleType_) <<
"\"", 1);
208 state->invPressureMassMatrix_ = getInvDiagonalOp(pressureMassMatrix_, scaleType_);
211 state->gamma_ = gamma_;
213 std::cout << Teuchos::describe(*(state->invPressureMassMatrix_), Teuchos::VERB_EXTREME) << std::endl;
217 if(state->B1tMpB1_ == Teuchos::null)
218 state->B1tMpB1_ = explicitMultiply(B1t, state->invPressureMassMatrix_, lpB1, state->B1tMpB1_);
222 LinearOp A11 = (rcp_dynamic_cast<const Thyra::DefaultAddedLinearOp<double> >(lpA11))->getOp(0);
223 state->A11p_ = explicitAdd(A11, scale(state->gamma_, state->B1tMpB1_), state->A11p_);
225 Teko_DEBUG_MSG(
"Computed A11p", 10);
227 if(state->B2tMpB2_ == Teuchos::null)
228 state->B2tMpB2_ = explicitMultiply(B2t, state->invPressureMassMatrix_, lpB2, state->B2tMpB2_);
229 LinearOp A22 = (rcp_dynamic_cast<const Thyra::DefaultAddedLinearOp<double> >(lpA22))->getOp(0);
230 state->A22p_ = explicitAdd(A22, scale(state->gamma_, state->B2tMpB2_), state->A22p_);
231 Teko_DEBUG_MSG(
"Computed A22p", 10);
235 if(state->B3tMpB3_ == Teuchos::null)
236 state->B3tMpB3_ = explicitMultiply(B3t, state->invPressureMassMatrix_, lpB3, state->B3tMpB3_);
237 LinearOp A33 = (rcp_dynamic_cast<const Thyra::DefaultAddedLinearOp<double> >(lpA33))->getOp(0);
238 state->A33p_ = explicitAdd(A33, scale(state->gamma_, state->B3tMpB3_), state->A33p_);
239 Teko_DEBUG_MSG(
"Computed A33p", 10);
243 if(state->isStabilized_)
245 if(state->S_ == Teuchos::null)
247 state->S_ = explicitAdd(scale(-1.0, lpC), scale(1.0/state->gamma_, pressureMassMatrix_), state->S_);
249 Teko_DEBUG_MSG(
"Computed S", 10);
252 state->setInitialized(
true);
256void InvModALStrategy::computeInverses(
const BlockedLinearOp & alOp,
257 ModALPrecondState *state)
const
260 TEUCHOS_ASSERT(dim == 2 || dim == 3);
262 Teko_DEBUG_SCOPE(
"InvModALStrategy::computeInverses", 10);
263 Teko_DEBUG_EXPR(Teuchos::Time invTimer(
""));
266 Teko_DEBUG_MSG(
"ModAL::computeInverses(): Building inv(A11)", 1);
267 Teko_DEBUG_EXPR(invTimer.start(
true));
269 InverseLinearOp invA11p = state->getInverse(
"invA11p");
270 if(invA11p == Teuchos::null)
273 state->addInverse(
"invA11p", invA11p);
280 Teko_DEBUG_EXPR(invTimer.stop());
281 Teko_DEBUG_MSG(
"ModAL::computeInverses GetInvA11 = " << invTimer.totalElapsedTime(), 1);
284 Teko_DEBUG_MSG(
"ModAL::computeInverses(): Building inv(A22)", 2);
285 Teko_DEBUG_EXPR(invTimer.start(
true));
287 InverseLinearOp invA22p = state->getInverse(
"invA22p");
288 if(invA22p == Teuchos::null)
291 state->addInverse(
"invA22p", invA22p);
298 Teko_DEBUG_EXPR(invTimer.stop());
299 Teko_DEBUG_MSG(
"ModAL::computeInverses(): GetInvA22 = " << invTimer.totalElapsedTime(), 2);
304 Teko_DEBUG_MSG(
"ModAL::computeInverses Building inv(A33)", 3);
305 Teko_DEBUG_EXPR(invTimer.start(
true));
307 InverseLinearOp invA33p = state->getInverse(
"invA33p");
308 if(invA33p == Teuchos::null)
311 state->addInverse(
"invA33p", invA33p);
318 Teko_DEBUG_EXPR(invTimer.stop());
319 Teko_DEBUG_MSG(
"ModAL::computeInverses GetInvA33 = " << invTimer.totalElapsedTime(), 3);
323 Teko_DEBUG_MSG(
"ModAL::computeInverses Building inv(S)", 4);
324 Teko_DEBUG_EXPR(invTimer.start(
true));
330 if(state->isStabilized_)
332 InverseLinearOp invS = state->getInverse(
"invS");
333 if(invS == Teuchos::null)
336 state->addInverse(
"invS", invS);
344 Teko_DEBUG_EXPR(invTimer.stop());
345 Teko_DEBUG_MSG(
"ModAL::computeInverses GetInvS = " << invTimer.totalElapsedTime(), 4);
InverseLinearOp buildInverse(const InverseFactory &factory, const LinearOp &A)
Build an inverse operator using a factory and a linear operator.
void rebuildInverse(const InverseFactory &factory, const LinearOp &A, InverseLinearOp &invA)
@ Diagonal
Specifies that just the diagonal is used.
int blockRowCount(const BlockedLinearOp &blo)
Get the row count in a block linear operator.
MultiVector getBlock(int i, const BlockedMultiVector &bmv)
Get the ith block from a BlockedMultiVector object.