Teko Version of the Day
Loading...
Searching...
No Matches
Teko_LSCSIMPLECStrategy.cpp
1/*
2// @HEADER
3//
4// ***********************************************************************
5//
6// Teko: A package for block and physics based preconditioning
7// Copyright 2010 Sandia Corporation
8//
9// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
10// the U.S. Government retains certain rights in this software.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric C. Cyr (eccyr@sandia.gov)
40//
41// ***********************************************************************
42//
43// @HEADER
44
45*/
46
47#include "Teko_LSCSIMPLECStrategy.hpp"
48
49#include "Thyra_DefaultDiagonalLinearOp.hpp"
50#include "Thyra_VectorStdOps.hpp"
51
52#include "Teuchos_Time.hpp"
53#include "Teuchos_TimeMonitor.hpp"
54
55// Teko includes
56#include "Teko_Utilities.hpp"
57#include "Teko_LSCPreconditionerFactory.hpp"
58
59using Teuchos::RCP;
60using Teuchos::rcp_dynamic_cast;
61using Teuchos::rcp_const_cast;
62
63namespace Teko {
64namespace NS {
65
67// LSCSIMPLECStrategy Implementation
69
70// constructors
72LSCSIMPLECStrategy::LSCSIMPLECStrategy()
73 : invFactoryF_(Teuchos::null), invFactoryS_(Teuchos::null)
74 , useFullLDU_(false), scaleType_(Diagonal)
75{ }
76
78
79void LSCSIMPLECStrategy::buildState(BlockedLinearOp & A,BlockPreconditionerState & state) const
80{
81 Teko_DEBUG_SCOPE("LSCSIMPLECStrategy::buildState",10);
82
83 LSCPrecondState * lscState = dynamic_cast<LSCPrecondState*>(&state);
84 TEUCHOS_ASSERT(lscState!=0);
85
86 // if neccessary save state information
87 if(not lscState->isInitialized()) {
88 Teko_DEBUG_EXPR(Teuchos::Time timer(""));
89
90 // construct operators
91 {
92 Teko_DEBUG_SCOPE("LSC-SIMPLEC::buildState constructing operators",1);
93 Teko_DEBUG_EXPR(timer.start(true));
94
95 initializeState(A,lscState);
96
97 Teko_DEBUG_EXPR(timer.stop());
98 Teko_DEBUG_MSG("LSC-SIMPLEC::buildState BuildOpsTime = " << timer.totalElapsedTime(),1);
99 }
100
101 // Build the inverses
102 {
103 Teko_DEBUG_SCOPE("LSC-SIMPLEC::buildState calculating inverses",1);
104 Teko_DEBUG_EXPR(timer.start(true));
105
106 computeInverses(A,lscState);
107
108 Teko_DEBUG_EXPR(timer.stop());
109 Teko_DEBUG_MSG("LSC-SIMPLEC::buildState BuildInvTime = " << timer.totalElapsedTime(),1);
110 }
111 }
112}
113
114// functions inherited from LSCStrategy
115LinearOp LSCSIMPLECStrategy::getInvBQBt(const BlockedLinearOp & /* A */,BlockPreconditionerState & state) const
116{
117 return state.getInverse("invBQBtmC");
118}
119
120LinearOp LSCSIMPLECStrategy::getInvBHBt(const BlockedLinearOp & /* A */,BlockPreconditionerState & state) const
121{
122 return state.getInverse("invBQBtmC").getConst();
123}
124
125LinearOp LSCSIMPLECStrategy::getInvF(const BlockedLinearOp & /* A */,BlockPreconditionerState & state) const
126{
127 return state.getInverse("invF");
128}
129
130LinearOp LSCSIMPLECStrategy::getInnerStabilization(const BlockedLinearOp & A,BlockPreconditionerState & state) const
131{
132 LSCPrecondState * lscState = dynamic_cast<LSCPrecondState*>(&state);
133 TEUCHOS_ASSERT(lscState!=0);
134 TEUCHOS_ASSERT(lscState->isInitialized())
135
136 const LinearOp C = getBlock(1,1,A);
137 return scale(-1.0,C);
138}
139
140
141LinearOp LSCSIMPLECStrategy::getInvMass(const BlockedLinearOp & /* A */,BlockPreconditionerState & state) const
142{
143 LSCPrecondState * lscState = dynamic_cast<LSCPrecondState*>(&state);
144 TEUCHOS_ASSERT(lscState!=0);
145 TEUCHOS_ASSERT(lscState->isInitialized())
146
147 return lscState->invMass_;
148}
149
150LinearOp LSCSIMPLECStrategy::getHScaling(const BlockedLinearOp & A,BlockPreconditionerState & state) const
151{
152 return getInvMass(A,state);
153}
154
156void LSCSIMPLECStrategy::initializeState(const BlockedLinearOp & A,LSCPrecondState * state) const
157{
158 Teko_DEBUG_SCOPE("LSCSIMPLECStrategy::initializeState",10);
159
160 const LinearOp F = getBlock(0,0,A);
161 const LinearOp Bt = getBlock(0,1,A);
162 const LinearOp B = getBlock(1,0,A);
163 const LinearOp C = getBlock(1,1,A);
164
165 bool isStabilized = (not isZeroOp(C));
166
167 state->invMass_ = getInvDiagonalOp(F,scaleType_);
168
169 // compute BQBt
170 state->BQBt_ = explicitMultiply(B,state->invMass_,Bt,state->BQBt_);
171 if(isStabilized) {
172 // now build B*Q*Bt-C
173 Teko::ModifiableLinearOp BQBtmC = state->getInverse("BQBtmC");
174 BQBtmC = explicitAdd(state->BQBt_,scale(-1.0,C),BQBtmC);
175 state->addInverse("BQBtmC",BQBtmC);
176 }
177 Teko_DEBUG_MSG("Computed BQBt",10);
178
179 state->setInitialized(true);
180}
181
187void LSCSIMPLECStrategy::computeInverses(const BlockedLinearOp & A,LSCPrecondState * state) const
188{
189 Teko_DEBUG_SCOPE("LSCSIMPLECStrategy::computeInverses",10);
190 Teko_DEBUG_EXPR(Teuchos::Time invTimer(""));
191
192 const LinearOp F = getBlock(0,0,A);
193
195
196 // (re)build the inverse of F
197 Teko_DEBUG_MSG("LSC-SIMPLEC::computeInverses Building inv(F)",1);
198 Teko_DEBUG_EXPR(invTimer.start(true));
199 InverseLinearOp invF = state->getInverse("invF");
200 if(invF==Teuchos::null) {
201 invF = buildInverse(*invFactoryF_,F);
202 state->addInverse("invF",invF);
203 } else {
204 rebuildInverse(*invFactoryF_,F,invF);
205 }
206 Teko_DEBUG_EXPR(invTimer.stop());
207 Teko_DEBUG_MSG("LSC-SIMPLEC::computeInverses GetInvF = " << invTimer.totalElapsedTime(),1);
208
210
211 // (re)build the inverse of BQBt
212 Teko_DEBUG_MSG("LSC-SIMPLEC::computeInverses Building inv(BQBtmC)",1);
213 Teko_DEBUG_EXPR(invTimer.start(true));
214 const LinearOp BQBt = state->getInverse("BQBtmC");
215 InverseLinearOp invBQBt = state->getInverse("invBQBtmC");
216 if(invBQBt==Teuchos::null) {
217 invBQBt = buildInverse(*invFactoryS_,BQBt);
218 state->addInverse("invBQBtmC",invBQBt);
219 } else {
220 rebuildInverse(*invFactoryS_,BQBt,invBQBt);
221 }
222 Teko_DEBUG_EXPR(invTimer.stop());
223 Teko_DEBUG_MSG("LSC-SIMPLEC::computeInverses GetInvBQBt = " << invTimer.totalElapsedTime(),1);
224}
225
227void LSCSIMPLECStrategy::initializeFromParameterList(const Teuchos::ParameterList & pl,const InverseLibrary & invLib)
228{
229 // get string specifying inverse
230 std::string invStr="", invVStr="", invPStr="";
231 bool useLDU = false;
232 scaleType_ = Diagonal;
233
234 // "parse" the parameter list
235 if(pl.isParameter("Inverse Type"))
236 invStr = pl.get<std::string>("Inverse Type");
237 if(pl.isParameter("Inverse Velocity Type"))
238 invVStr = pl.get<std::string>("Inverse Velocity Type");
239 if(pl.isParameter("Inverse Pressure Type"))
240 invPStr = pl.get<std::string>("Inverse Pressure Type");
241 if(pl.isParameter("Use LDU"))
242 useLDU = pl.get<bool>("Use LDU");
243 if(pl.isParameter("Scaling Type")) {
244 scaleType_ = getDiagonalType(pl.get<std::string>("Scaling Type"));
245 TEUCHOS_TEST_FOR_EXCEPT(scaleType_==NotDiag);
246 }
247
248 Teko_DEBUG_MSG_BEGIN(0)
249 DEBUG_STREAM << "LSC Inverse Strategy Parameters: " << std::endl;
250 DEBUG_STREAM << " inv type = \"" << invStr << "\"" << std::endl;
251 DEBUG_STREAM << " inv v type = \"" << invVStr << "\"" << std::endl;
252 DEBUG_STREAM << " inv p type = \"" << invPStr << "\"" << std::endl;
253 DEBUG_STREAM << " use ldu = " << useLDU << std::endl;
254 DEBUG_STREAM << " scale type = " << getDiagonalName(scaleType_) << std::endl;
255 DEBUG_STREAM << "LSC Inverse Strategy Parameter list: " << std::endl;
256 pl.print(DEBUG_STREAM);
257 Teko_DEBUG_MSG_END()
258
259 // set defaults as needed
260 if(invStr=="") invStr = "Amesos";
261 if(invVStr=="") invVStr = invStr;
262 if(invPStr=="") invPStr = invStr;
263
264 // build velocity inverse factory
265 invFactoryF_ = invLib.getInverseFactory(invVStr);
266 invFactoryS_ = invFactoryF_; // by default these are the same
267 if(invVStr!=invPStr) // if different, build pressure inverse factory
268 invFactoryS_ = invLib.getInverseFactory(invPStr);
269
270 // set other parameters
271 setUseFullLDU(useLDU);
272}
273
274} // end namespace NS
275} // 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)
@ NotDiag
For user convenience, if Teko recieves this value, exceptions will be thrown.
@ Diagonal
Specifies that just the diagonal is used.
MultiVector getBlock(int i, const BlockedMultiVector &bmv)
Get the ith block from a BlockedMultiVector object.
An implementation of a state object for block preconditioners.
Preconditioner state for the LSC factory.