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"
19 using Teuchos::rcp_const_cast;
20 using Teuchos::rcp_dynamic_cast;
27 InvModALStrategy::InvModALStrategy()
28 : invFactoryA_(Teuchos::null),
29 invFactoryS_(Teuchos::null),
30 pressureMassMatrix_(Teuchos::null),
36 InvModALStrategy::InvModALStrategy(
const Teuchos::RCP<InverseFactory> &factory)
37 : invFactoryA_(factory),
38 invFactoryS_(factory),
39 pressureMassMatrix_(Teuchos::null),
45 InvModALStrategy::InvModALStrategy(
const Teuchos::RCP<InverseFactory> &invFactoryA,
46 const Teuchos::RCP<InverseFactory> &invFactoryS)
47 : invFactoryA_(invFactoryA),
48 invFactoryS_(invFactoryS),
49 pressureMassMatrix_(Teuchos::null),
55 InvModALStrategy::InvModALStrategy(
const Teuchos::RCP<InverseFactory> &invFactory,
56 LinearOp &pressureMassMatrix)
57 : invFactoryA_(invFactory),
58 invFactoryS_(invFactory),
59 pressureMassMatrix_(pressureMassMatrix),
65 InvModALStrategy::InvModALStrategy(
const Teuchos::RCP<InverseFactory> &invFactoryA,
66 const Teuchos::RCP<InverseFactory> &invFactoryS,
67 LinearOp &pressureMassMatrix)
68 : invFactoryA_(invFactoryA),
69 invFactoryS_(invFactoryS),
70 pressureMassMatrix_(pressureMassMatrix),
76 LinearOp InvModALStrategy::getInvA11p(BlockPreconditionerState &state)
const {
77 return state.getInverse(
"invA11p");
80 LinearOp InvModALStrategy::getInvA22p(BlockPreconditionerState &state)
const {
81 return state.getInverse(
"invA22p");
84 LinearOp InvModALStrategy::getInvA33p(BlockPreconditionerState &state)
const {
85 return state.getInverse(
"invA33p");
88 LinearOp InvModALStrategy::getInvS(BlockPreconditionerState &state)
const {
89 return state.getInverse(
"invS");
93 void InvModALStrategy::setPressureMassMatrix(
const LinearOp &pressureMassMatrix) {
94 pressureMassMatrix_ = pressureMassMatrix;
98 void InvModALStrategy::setGamma(
double gamma) {
99 TEUCHOS_ASSERT(gamma > 0.0);
103 void InvModALStrategy::buildState(
const BlockedLinearOp &alOp,
104 BlockPreconditionerState &state)
const {
105 Teko_DEBUG_SCOPE(
"InvModALStrategy::buildState", 10);
107 ModALPrecondState *modALState =
dynamic_cast<ModALPrecondState *
>(&state);
108 TEUCHOS_ASSERT(modALState != NULL);
111 if (not modALState->isInitialized()) {
112 Teko_DEBUG_EXPR(Teuchos::Time timer(
""));
116 Teko_DEBUG_SCOPE(
"ModAL::buildState: Initializing state object", 1);
117 Teko_DEBUG_EXPR(timer.start(
true));
119 initializeState(alOp, modALState);
121 Teko_DEBUG_EXPR(timer.stop());
122 Teko_DEBUG_MSG(
"ModAL::buildState: BuildOpsTime = " << timer.totalElapsedTime(), 1);
127 Teko_DEBUG_SCOPE(
"ModAL::buildState: Computing inverses", 1);
128 Teko_DEBUG_EXPR(timer.start(
true));
130 computeInverses(alOp, modALState);
132 Teko_DEBUG_EXPR(timer.stop());
133 Teko_DEBUG_MSG(
"ModAL::buildState: BuildInvTime = " << timer.totalElapsedTime(), 1);
139 void InvModALStrategy::initializeState(
const BlockedLinearOp &alOp,
140 ModALPrecondState *state)
const {
141 Teko_DEBUG_SCOPE(
"InvModALStrategy::initializeState", 10);
144 int dim = blockRowCount(alOp) - 1;
145 TEUCHOS_ASSERT(dim == 2 || dim == 3);
147 LinearOp lpA11 = getBlock(0, 0, alOp);
148 LinearOp lpA22 = getBlock(1, 1, alOp);
149 LinearOp lpA33, lpB1, lpB2, lpB3, lpB1t, lpB2t, lpB3t, lpC;
153 lpB1 = getBlock(2, 0, alOp);
154 lpB2 = getBlock(2, 1, alOp);
155 lpB1t = getBlock(0, 2, alOp);
156 lpB2t = getBlock(1, 2, alOp);
157 lpC = getBlock(2, 2, alOp);
161 lpA33 = getBlock(2, 2, alOp);
162 lpB1 = getBlock(3, 0, alOp);
163 lpB2 = getBlock(3, 1, alOp);
164 lpB3 = getBlock(3, 2, alOp);
165 lpB1t = getBlock(0, 3, alOp);
166 lpB2t = getBlock(1, 3, alOp);
167 lpB3t = getBlock(2, 3, alOp);
168 lpC = getBlock(3, 3, alOp);
173 LinearOp B1t = (rcp_dynamic_cast<
const Thyra::DefaultAddedLinearOp<double> >(lpB1t))->getOp(0);
174 LinearOp B2t = (rcp_dynamic_cast<
const Thyra::DefaultAddedLinearOp<double> >(lpB2t))->getOp(0);
177 B3t = (rcp_dynamic_cast<
const Thyra::DefaultAddedLinearOp<double> >(lpB3t))->getOp(0);
182 state->isStabilized_ = (not isZeroOp(lpC));
186 state->pressureMassMatrix_ = pressureMassMatrix_;
188 if (state->pressureMassMatrix_ == Teuchos::null) {
189 Teko_DEBUG_MSG(
"InvModALStrategy::initializeState(): Build identity type \""
190 << getDiagonalName(scaleType_) <<
"\"",
192 state->invPressureMassMatrix_ = Thyra::identity<double>(lpB1->range());
196 else if (state->invPressureMassMatrix_ == Teuchos::null) {
197 Teko_DEBUG_MSG(
"ModAL::initializeState(): Build Scaling <mass> type \""
198 << getDiagonalName(scaleType_) <<
"\"",
200 state->invPressureMassMatrix_ = getInvDiagonalOp(pressureMassMatrix_, scaleType_);
203 state->gamma_ = gamma_;
205 std::cout << Teuchos::describe(*(state->invPressureMassMatrix_), Teuchos::VERB_EXTREME)
210 if (state->B1tMpB1_ == Teuchos::null)
211 state->B1tMpB1_ = explicitMultiply(B1t, state->invPressureMassMatrix_, lpB1, state->B1tMpB1_);
215 LinearOp A11 = (rcp_dynamic_cast<
const Thyra::DefaultAddedLinearOp<double> >(lpA11))->getOp(0);
216 state->A11p_ = explicitAdd(A11, scale(state->gamma_, state->B1tMpB1_), state->A11p_);
218 Teko_DEBUG_MSG(
"Computed A11p", 10);
220 if (state->B2tMpB2_ == Teuchos::null)
221 state->B2tMpB2_ = explicitMultiply(B2t, state->invPressureMassMatrix_, lpB2, state->B2tMpB2_);
222 LinearOp A22 = (rcp_dynamic_cast<
const Thyra::DefaultAddedLinearOp<double> >(lpA22))->getOp(0);
223 state->A22p_ = explicitAdd(A22, scale(state->gamma_, state->B2tMpB2_), state->A22p_);
224 Teko_DEBUG_MSG(
"Computed A22p", 10);
227 if (state->B3tMpB3_ == Teuchos::null)
228 state->B3tMpB3_ = explicitMultiply(B3t, state->invPressureMassMatrix_, lpB3, state->B3tMpB3_);
229 LinearOp A33 = (rcp_dynamic_cast<
const Thyra::DefaultAddedLinearOp<double> >(lpA33))->getOp(0);
230 state->A33p_ = explicitAdd(A33, scale(state->gamma_, state->B3tMpB3_), state->A33p_);
231 Teko_DEBUG_MSG(
"Computed A33p", 10);
235 if (state->isStabilized_) {
236 if (state->S_ == Teuchos::null) {
238 explicitAdd(scale(-1.0, lpC), scale(1.0 / state->gamma_, pressureMassMatrix_), state->S_);
240 Teko_DEBUG_MSG(
"Computed S", 10);
243 state->setInitialized(
true);
247 void InvModALStrategy::computeInverses(
const BlockedLinearOp &alOp,
248 ModALPrecondState *state)
const {
249 int dim = blockRowCount(alOp) - 1;
250 TEUCHOS_ASSERT(dim == 2 || dim == 3);
252 Teko_DEBUG_SCOPE(
"InvModALStrategy::computeInverses", 10);
253 Teko_DEBUG_EXPR(Teuchos::Time invTimer(
""));
256 Teko_DEBUG_MSG(
"ModAL::computeInverses(): Building inv(A11)", 1);
257 Teko_DEBUG_EXPR(invTimer.start(
true));
259 InverseLinearOp invA11p = state->getInverse(
"invA11p");
260 if (invA11p == Teuchos::null) {
262 state->addInverse(
"invA11p", invA11p);
267 Teko_DEBUG_EXPR(invTimer.stop());
268 Teko_DEBUG_MSG(
"ModAL::computeInverses GetInvA11 = " << invTimer.totalElapsedTime(), 1);
271 Teko_DEBUG_MSG(
"ModAL::computeInverses(): Building inv(A22)", 2);
272 Teko_DEBUG_EXPR(invTimer.start(
true));
274 InverseLinearOp invA22p = state->getInverse(
"invA22p");
275 if (invA22p == Teuchos::null) {
277 state->addInverse(
"invA22p", invA22p);
282 Teko_DEBUG_EXPR(invTimer.stop());
283 Teko_DEBUG_MSG(
"ModAL::computeInverses(): GetInvA22 = " << invTimer.totalElapsedTime(), 2);
287 Teko_DEBUG_MSG(
"ModAL::computeInverses Building inv(A33)", 3);
288 Teko_DEBUG_EXPR(invTimer.start(
true));
290 InverseLinearOp invA33p = state->getInverse(
"invA33p");
291 if (invA33p == Teuchos::null) {
293 state->addInverse(
"invA33p", invA33p);
298 Teko_DEBUG_EXPR(invTimer.stop());
299 Teko_DEBUG_MSG(
"ModAL::computeInverses GetInvA33 = " << invTimer.totalElapsedTime(), 3);
303 Teko_DEBUG_MSG(
"ModAL::computeInverses Building inv(S)", 4);
304 Teko_DEBUG_EXPR(invTimer.start(
true));
310 if (state->isStabilized_) {
311 InverseLinearOp invS = state->getInverse(
"invS");
312 if (invS == Teuchos::null) {
314 state->addInverse(
"invS", invS);
320 Teko_DEBUG_EXPR(invTimer.stop());
321 Teko_DEBUG_MSG(
"ModAL::computeInverses GetInvS = " << invTimer.totalElapsedTime(), 4);
void rebuildInverse(const InverseFactory &factory, const LinearOp &A, InverseLinearOp &invA)
InverseLinearOp buildInverse(const InverseFactory &factory, const LinearOp &A)
Build an inverse operator using a factory and a linear operator.