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