Teko  Version of the Day
 All Classes Files Functions Variables Pages
Teko_InvModALStrategy.cpp
1 /*
2  * Author: Zhen Wang
3  * Email: wangz@ornl.gov
4  * zhen.wang@alum.emory.edu
5  */
6 
7 #include "Thyra_DefaultDiagonalLinearOp.hpp"
8 #include "Thyra_DefaultAddedLinearOp_def.hpp"
9 #include "Thyra_DefaultIdentityLinearOp_decl.hpp"
10 
11 #include "Teuchos_Time.hpp"
12 
13 #include "Teko_Utilities.hpp"
14 
15 #include "Teko_InvModALStrategy.hpp"
16 #include "Teko_ModALPreconditionerFactory.hpp"
17 
18 using Teuchos::RCP;
19 using Teuchos::rcp_const_cast;
20 using Teuchos::rcp_dynamic_cast;
21 
22 namespace Teko {
23 
24 namespace NS {
25 
26 // Empty constructor.
27 InvModALStrategy::InvModALStrategy()
28  : invFactoryA_(Teuchos::null),
29  invFactoryS_(Teuchos::null),
30  pressureMassMatrix_(Teuchos::null),
31  gamma_(0.05),
32  scaleType_(Diagonal),
33  isSymmetric_(true) {}
34 
35 // If there is only one InverseFactory, use it for all solves.
36 InvModALStrategy::InvModALStrategy(const Teuchos::RCP<InverseFactory> &factory)
37  : invFactoryA_(factory),
38  invFactoryS_(factory),
39  pressureMassMatrix_(Teuchos::null),
40  gamma_(0.05),
41  scaleType_(Diagonal),
42  isSymmetric_(true) {}
43 
44 // If there are two InverseFactory's...
45 InvModALStrategy::InvModALStrategy(const Teuchos::RCP<InverseFactory> &invFactoryA,
46  const Teuchos::RCP<InverseFactory> &invFactoryS)
47  : invFactoryA_(invFactoryA),
48  invFactoryS_(invFactoryS),
49  pressureMassMatrix_(Teuchos::null),
50  gamma_(0.05),
51  scaleType_(Diagonal),
52  isSymmetric_(true) {}
53 
54 // If there are two InverseFactory's...
55 InvModALStrategy::InvModALStrategy(const Teuchos::RCP<InverseFactory> &invFactory,
56  LinearOp &pressureMassMatrix)
57  : invFactoryA_(invFactory),
58  invFactoryS_(invFactory),
59  pressureMassMatrix_(pressureMassMatrix),
60  gamma_(0.05),
61  scaleType_(Diagonal),
62  isSymmetric_(true) {}
63 
64 // If there are two InverseFactory's...
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),
71  gamma_(0.05),
72  scaleType_(Diagonal),
73  isSymmetric_(true) {}
74 
75 // Return "inverses".
76 LinearOp InvModALStrategy::getInvA11p(BlockPreconditionerState &state) const {
77  return state.getInverse("invA11p");
78 }
79 
80 LinearOp InvModALStrategy::getInvA22p(BlockPreconditionerState &state) const {
81  return state.getInverse("invA22p");
82 }
83 
84 LinearOp InvModALStrategy::getInvA33p(BlockPreconditionerState &state) const {
85  return state.getInverse("invA33p");
86 }
87 
88 LinearOp InvModALStrategy::getInvS(BlockPreconditionerState &state) const {
89  return state.getInverse("invS");
90 }
91 
92 // Set pressure mass matrix.
93 void InvModALStrategy::setPressureMassMatrix(const LinearOp &pressureMassMatrix) {
94  pressureMassMatrix_ = pressureMassMatrix;
95 }
96 
97 // Set gamma.
98 void InvModALStrategy::setGamma(double gamma) {
99  TEUCHOS_ASSERT(gamma > 0.0);
100  gamma_ = gamma;
101 }
102 
103 void InvModALStrategy::buildState(const BlockedLinearOp &alOp,
104  BlockPreconditionerState &state) const {
105  Teko_DEBUG_SCOPE("InvModALStrategy::buildState", 10);
106 
107  ModALPrecondState *modALState = dynamic_cast<ModALPrecondState *>(&state);
108  TEUCHOS_ASSERT(modALState != NULL);
109 
110  // if necessary save state information
111  if (not modALState->isInitialized()) {
112  Teko_DEBUG_EXPR(Teuchos::Time timer(""));
113 
114  {
115  // construct operators
116  Teko_DEBUG_SCOPE("ModAL::buildState: Initializing state object", 1);
117  Teko_DEBUG_EXPR(timer.start(true));
118 
119  initializeState(alOp, modALState);
120 
121  Teko_DEBUG_EXPR(timer.stop());
122  Teko_DEBUG_MSG("ModAL::buildState: BuildOpsTime = " << timer.totalElapsedTime(), 1);
123  }
124 
125  {
126  // Build the inverses
127  Teko_DEBUG_SCOPE("ModAL::buildState: Computing inverses", 1);
128  Teko_DEBUG_EXPR(timer.start(true));
129 
130  computeInverses(alOp, modALState);
131 
132  Teko_DEBUG_EXPR(timer.stop());
133  Teko_DEBUG_MSG("ModAL::buildState: BuildInvTime = " << timer.totalElapsedTime(), 1);
134  }
135  }
136 }
137 
138 // Initialize the state object using the ALOperator.
139 void InvModALStrategy::initializeState(const BlockedLinearOp &alOp,
140  ModALPrecondState *state) const {
141  Teko_DEBUG_SCOPE("InvModALStrategy::initializeState", 10);
142 
143  // Extract sub-matrices from blocked linear operator.
144  int dim = blockRowCount(alOp) - 1;
145  TEUCHOS_ASSERT(dim == 2 || dim == 3);
146 
147  LinearOp lpA11 = getBlock(0, 0, alOp);
148  LinearOp lpA22 = getBlock(1, 1, alOp);
149  LinearOp lpA33, lpB1, lpB2, lpB3, lpB1t, lpB2t, lpB3t, lpC;
150 
151  // 2D problem.
152  if (dim == 2) {
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);
158  }
159  // 3D problem.
160  else if (dim == 3) {
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);
169  }
170 
171  // For problems using stabilized finite elements,
172  // lpB1t, lpB2t and lpB3t are added linear operators. Extract original operators.
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);
175  LinearOp B3t;
176  if (dim == 3) {
177  B3t = (rcp_dynamic_cast<const Thyra::DefaultAddedLinearOp<double> >(lpB3t))->getOp(0);
178  }
179 
180  // std::cout << Teuchos::describe(*lpC, Teuchos::VERB_EXTREME) << std::endl;
181  // Check whether the finite elements are stable or not.
182  state->isStabilized_ = (not isZeroOp(lpC));
183  // state->isStabilized_ = false;
184  // std::cout << state->isStabilized_ << std::endl;
185 
186  state->pressureMassMatrix_ = pressureMassMatrix_;
187  // If pressure mass matrix is not set, use identity.
188  if (state->pressureMassMatrix_ == Teuchos::null) {
189  Teko_DEBUG_MSG("InvModALStrategy::initializeState(): Build identity type \""
190  << getDiagonalName(scaleType_) << "\"",
191  1);
192  state->invPressureMassMatrix_ = Thyra::identity<double>(lpB1->range());
193  }
194  // If the inverse of the pressure mass matrix is not set,
195  // build it from the pressure mass matrix.
196  else if (state->invPressureMassMatrix_ == Teuchos::null) {
197  Teko_DEBUG_MSG("ModAL::initializeState(): Build Scaling <mass> type \""
198  << getDiagonalName(scaleType_) << "\"",
199  1);
200  state->invPressureMassMatrix_ = getInvDiagonalOp(pressureMassMatrix_, scaleType_);
201  }
202  // Else "invPressureMassMatrix_" should be set and there is no reason to rebuild it
203  state->gamma_ = gamma_;
204  // S_ = scale(1.0/gamma_, pressureMassMatrix_);
205  std::cout << Teuchos::describe(*(state->invPressureMassMatrix_), Teuchos::VERB_EXTREME)
206  << std::endl;
207 
208  // Build state variables: B_1^T*W^{-1}*B_1, A11p, etc.
209  // B_1^T*W^{-1}*B_1 may not change so save it in the state.
210  if (state->B1tMpB1_ == Teuchos::null)
211  state->B1tMpB1_ = explicitMultiply(B1t, state->invPressureMassMatrix_, lpB1, state->B1tMpB1_);
212  // Get the(1,1) block of the non-augmented matrix.
213  // Recall alOp is augmented. So lpA11 = A11 + gamma B_1^T W^{-1} B_1.
214  // Cast lpA11 as an added linear operator and get the first item.
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_);
217  // std::cout << Teuchos::describe(*(state->B1tMpB1_), Teuchos::VERB_EXTREME) << std::endl;
218  Teko_DEBUG_MSG("Computed A11p", 10);
219 
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);
225 
226  if (dim == 3) {
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);
232  }
233 
234  // Inverse the Schur complement.
235  if (state->isStabilized_) {
236  if (state->S_ == Teuchos::null) {
237  state->S_ =
238  explicitAdd(scale(-1.0, lpC), scale(1.0 / state->gamma_, pressureMassMatrix_), state->S_);
239  }
240  Teko_DEBUG_MSG("Computed S", 10);
241  }
242 
243  state->setInitialized(true);
244 }
245 
246 // Compute inverses.
247 void InvModALStrategy::computeInverses(const BlockedLinearOp &alOp,
248  ModALPrecondState *state) const {
249  int dim = blockRowCount(alOp) - 1;
250  TEUCHOS_ASSERT(dim == 2 || dim == 3);
251 
252  Teko_DEBUG_SCOPE("InvModALStrategy::computeInverses", 10);
253  Teko_DEBUG_EXPR(Teuchos::Time invTimer(""));
254 
255  //(re)build the inverse of A11
256  Teko_DEBUG_MSG("ModAL::computeInverses(): Building inv(A11)", 1);
257  Teko_DEBUG_EXPR(invTimer.start(true));
258 
259  InverseLinearOp invA11p = state->getInverse("invA11p");
260  if (invA11p == Teuchos::null) {
261  invA11p = buildInverse(*invFactoryA_, state->A11p_);
262  state->addInverse("invA11p", invA11p);
263  } else {
264  rebuildInverse(*invFactoryA_, state->A11p_, invA11p);
265  }
266 
267  Teko_DEBUG_EXPR(invTimer.stop());
268  Teko_DEBUG_MSG("ModAL::computeInverses GetInvA11 = " << invTimer.totalElapsedTime(), 1);
269 
270  //(re)build the inverse of A22
271  Teko_DEBUG_MSG("ModAL::computeInverses(): Building inv(A22)", 2);
272  Teko_DEBUG_EXPR(invTimer.start(true));
273 
274  InverseLinearOp invA22p = state->getInverse("invA22p");
275  if (invA22p == Teuchos::null) {
276  invA22p = buildInverse(*invFactoryA_, state->A22p_);
277  state->addInverse("invA22p", invA22p);
278  } else {
279  rebuildInverse(*invFactoryA_, state->A22p_, invA22p);
280  }
281 
282  Teko_DEBUG_EXPR(invTimer.stop());
283  Teko_DEBUG_MSG("ModAL::computeInverses(): GetInvA22 = " << invTimer.totalElapsedTime(), 2);
284 
285  //(re)build the inverse of A33
286  if (dim == 3) {
287  Teko_DEBUG_MSG("ModAL::computeInverses Building inv(A33)", 3);
288  Teko_DEBUG_EXPR(invTimer.start(true));
289 
290  InverseLinearOp invA33p = state->getInverse("invA33p");
291  if (invA33p == Teuchos::null) {
292  invA33p = buildInverse(*invFactoryA_, state->A33p_);
293  state->addInverse("invA33p", invA33p);
294  } else {
295  rebuildInverse(*invFactoryA_, state->A33p_, invA33p);
296  }
297 
298  Teko_DEBUG_EXPR(invTimer.stop());
299  Teko_DEBUG_MSG("ModAL::computeInverses GetInvA33 = " << invTimer.totalElapsedTime(), 3);
300  }
301 
302  //(re)build the inverse of S
303  Teko_DEBUG_MSG("ModAL::computeInverses Building inv(S)", 4);
304  Teko_DEBUG_EXPR(invTimer.start(true));
305 
306  // There are two ways to "invert" S.
307  // The following method construct invS by InverseFactory invFactoryS_.
308  // The other way is to use diagonal approximation,
309  // which is done in ModALPreconditionerFactory.cpp.
310  if (state->isStabilized_) {
311  InverseLinearOp invS = state->getInverse("invS");
312  if (invS == Teuchos::null) {
313  invS = buildInverse(*invFactoryS_, state->S_);
314  state->addInverse("invS", invS);
315  } else {
316  rebuildInverse(*invFactoryS_, state->S_, invS);
317  }
318  }
319 
320  Teko_DEBUG_EXPR(invTimer.stop());
321  Teko_DEBUG_MSG("ModAL::computeInverses GetInvS = " << invTimer.totalElapsedTime(), 4);
322 }
323 
324 } // end namespace NS
325 
326 } // end namespace Teko
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.