Teko Version of the Day
Loading...
Searching...
No Matches
Teko_PresLaplaceLSCStrategy.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 "NS/Teko_PresLaplaceLSCStrategy.hpp"
48
49#include "Thyra_DefaultDiagonalLinearOp.hpp"
50
51#include "Teuchos_Time.hpp"
52#include "Teuchos_TimeMonitor.hpp"
53
54// Teko includes
55#include "Teko_Utilities.hpp"
56#include "NS/Teko_LSCPreconditionerFactory.hpp"
57
58using Teuchos::RCP;
59using Teuchos::rcp_dynamic_cast;
60using Teuchos::rcp_const_cast;
61
62namespace Teko {
63namespace NS {
64
66// PresLaplaceLSCStrategy Implementation
68
69// constructors
71PresLaplaceLSCStrategy::PresLaplaceLSCStrategy()
72 : invFactoryV_(Teuchos::null), invFactoryP_(Teuchos::null)
73 , eigSolveParam_(5), useFullLDU_(false), useMass_(false), scaleType_(AbsRowSum)
74{ }
75
76PresLaplaceLSCStrategy::PresLaplaceLSCStrategy(const Teuchos::RCP<InverseFactory> & factory)
77 : invFactoryV_(factory), invFactoryP_(factory)
78 , eigSolveParam_(5), useFullLDU_(false), useMass_(false), scaleType_(AbsRowSum)
79{ }
80
81PresLaplaceLSCStrategy::PresLaplaceLSCStrategy(const Teuchos::RCP<InverseFactory> & invFactF,
82 const Teuchos::RCP<InverseFactory> & invFactS)
83 : invFactoryV_(invFactF), invFactoryP_(invFactS)
84 , eigSolveParam_(5), useFullLDU_(false), useMass_(false), scaleType_(AbsRowSum)
85{ }
86
88
89void PresLaplaceLSCStrategy::buildState(BlockedLinearOp & A,BlockPreconditionerState & state) const
90{
91 Teko_DEBUG_SCOPE("PresLaplaceLSCStrategy::buildState",10);
92
93 LSCPrecondState * lscState = dynamic_cast<LSCPrecondState*>(&state);
94 TEUCHOS_ASSERT(lscState!=0);
95
96 // if neccessary save state information
97 if(not lscState->isInitialized()) {
98 Teko_DEBUG_EXPR(Teuchos::Time timer(""));
99
100 // construct operators
101 {
102 Teko_DEBUG_SCOPE("PL-LSC::buildState constructing operators",1);
103 Teko_DEBUG_EXPR(timer.start(true));
104
105 initializeState(A,lscState);
106
107 Teko_DEBUG_EXPR(timer.stop());
108 Teko_DEBUG_MSG("PL-LSC::buildState BuildOpsTime = " << timer.totalElapsedTime(),1);
109 }
110
111 // Build the inverses
112 {
113 Teko_DEBUG_SCOPE("PL-LSC::buildState calculating inverses",1);
114 Teko_DEBUG_EXPR(timer.start(true));
115
116 computeInverses(A,lscState);
117
118 Teko_DEBUG_EXPR(timer.stop());
119 Teko_DEBUG_MSG("PL-LSC::buildState BuildInvTime = " << timer.totalElapsedTime(),1);
120 }
121 }
122}
123
124// functions inherited from LSCStrategy
125LinearOp PresLaplaceLSCStrategy::getInvBQBt(const BlockedLinearOp & /* A */,BlockPreconditionerState & state) const
126{
127 return state.getModifiableOp("invPresLap");
128}
129
130LinearOp PresLaplaceLSCStrategy::getInvBHBt(const BlockedLinearOp & /* A */,BlockPreconditionerState & state) const
131{
132 return state.getModifiableOp("invPresLap");
133}
134
135LinearOp PresLaplaceLSCStrategy::getInvF(const BlockedLinearOp & /* A */,BlockPreconditionerState & state) const
136{
137 return state.getModifiableOp("invF");
138}
139
140// LinearOp PresLaplaceLSCStrategy::getInvAlphaD(const BlockedLinearOp & A,BlockPreconditionerState & state) const
141LinearOp PresLaplaceLSCStrategy::getOuterStabilization(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->aiD_;
148}
149
150LinearOp PresLaplaceLSCStrategy::getInvMass(const BlockedLinearOp & /* A */,BlockPreconditionerState & state) const
151{
152 LSCPrecondState * lscState = dynamic_cast<LSCPrecondState*>(&state);
153 TEUCHOS_ASSERT(lscState!=0);
154 TEUCHOS_ASSERT(lscState->isInitialized())
155
156 return lscState->invMass_;
157}
158
159LinearOp PresLaplaceLSCStrategy::getHScaling(const BlockedLinearOp & A,BlockPreconditionerState & state) const
160{
161 return getInvMass(A,state);
162}
163
165void PresLaplaceLSCStrategy::initializeState(const BlockedLinearOp & A,LSCPrecondState * state) const
166{
167 Teko_DEBUG_SCOPE("PresLaplaceLSCStrategy::initializeState",10);
168
169 std::string velMassStr = getVelocityMassString();
170
171 const LinearOp F = getBlock(0,0,A);
172 const LinearOp Bt = getBlock(0,1,A);
173 const LinearOp B = getBlock(1,0,A);
174 const LinearOp C = getBlock(1,1,A);
175
176 LinearOp D = B;
177 LinearOp G = Bt;
178
179 bool isStabilized = (not isZeroOp(C));
180
181 // grab operators from state object
182 LinearOp massMatrix = state->getLinearOp(velMassStr);
183
184 // The logic follows like this
185 // if there is no mass matrix available --> build from F
186 // if there is a mass matrix and the inverse hasn't yet been built
187 // --> build from the mass matrix
188 // otherwise, there is already an invMass_ matrix that is appropriate
189 // --> use that one
190 if(massMatrix==Teuchos::null) {
191 Teko_DEBUG_MSG("PL-LSC::initializeState Build Scaling <F> type \""
192 << getDiagonalName(scaleType_) << "\"" ,1);
193 state->invMass_ = getInvDiagonalOp(F,scaleType_);
194 }
195 else if(state->invMass_==Teuchos::null) {
196 Teko_DEBUG_MSG("PL-LSC::initializeState Build Scaling <mass> type \""
197 << getDiagonalName(scaleType_) << "\"" ,1);
198 state->invMass_ = getInvDiagonalOp(massMatrix,scaleType_);
199 }
200 // else "invMass_" should be set and there is no reason to rebuild it
201
202 // if this is a stable discretization...we are done!
203 if(not isStabilized) {
204 state->aiD_ = Teuchos::null;
205
206 state->setInitialized(true);
207
208 return;
209 }
210
211 // compute alpha scaled inv(D): EHSST2007 Eq. 4.29
212 // construct B_idF_Bt and save it for refilling later: This could reuse BQBt graph
213 LinearOp invDiagF = getInvDiagonalOp(F);
214 Teko::ModifiableLinearOp & modB_idF_Bt = state->getModifiableOp("BidFBt");
215 modB_idF_Bt = explicitMultiply(B,invDiagF,Bt,modB_idF_Bt);
216 const LinearOp B_idF_Bt = modB_idF_Bt;
217
218 MultiVector vec_D = getDiagonal(B_idF_Bt); // this memory could be reused
219 update(-1.0,getDiagonal(C),1.0,vec_D); // vec_D = diag(B*inv(diag(F))*Bt)-diag(C)
220 const LinearOp invD = buildInvDiagonal(vec_D,"inv(D)");
221
222 Teko_DEBUG_MSG("Calculating alpha",10);
223 const LinearOp BidFBtidD = multiply<double>(B_idF_Bt,invD);
224 double num = std::fabs(Teko::computeSpectralRad(BidFBtidD,5e-2,false,eigSolveParam_));
225 Teko_DEBUG_MSG("Calculated alpha",10);
226 state->alpha_ = 1.0/num;
227 state->aiD_ = Thyra::scale(state->alpha_,invD);
228
229 Teko_DEBUG_MSG_BEGIN(5)
230 DEBUG_STREAM << "PL-LSC Alpha Parameter = " << state->alpha_ << std::endl;
231 Teko_DEBUG_MSG_END()
232
233 state->setInitialized(true);
234}
235
241void PresLaplaceLSCStrategy::computeInverses(const BlockedLinearOp & A,LSCPrecondState * state) const
242{
243 Teko_DEBUG_SCOPE("PresLaplaceLSCStrategy::computeInverses",10);
244 Teko_DEBUG_EXPR(Teuchos::Time invTimer(""));
245
246 std::string presLapStr = getPressureLaplaceString();
247
248 const LinearOp F = getBlock(0,0,A);
249 const LinearOp presLap = state->getLinearOp(presLapStr);
250
252
253 // (re)build the inverse of F
254 Teko_DEBUG_MSG("PL-LSC::computeInverses Building inv(F)",1);
255 Teko_DEBUG_EXPR(invTimer.start(true));
256 ModifiableLinearOp & invF = state->getModifiableOp("invF");
257 if(invF==Teuchos::null) {
258 invF = buildInverse(*invFactoryV_,F);
259 } else {
260 rebuildInverse(*invFactoryV_,F,invF);
261 }
262 Teko_DEBUG_EXPR(invTimer.stop());
263 Teko_DEBUG_MSG("PL-LSC::computeInverses GetInvF = " << invTimer.totalElapsedTime(),1);
264
266
267 // (re)build the inverse of P
268 Teko_DEBUG_MSG("PL-LSC::computeInverses Building inv(PresLap)",1);
269 Teko_DEBUG_EXPR(invTimer.start(true));
270 ModifiableLinearOp & invPresLap = state->getModifiableOp("invPresLap");
271 if(invPresLap==Teuchos::null) {
272 invPresLap = buildInverse(*invFactoryP_,presLap);
273 } else {
274 // not need because the pressure laplacian never changes
275 // rebuildInverse(*invFactoryP_,presLap,invPresLap);
276 }
277 Teko_DEBUG_EXPR(invTimer.stop());
278 Teko_DEBUG_MSG("PL-LSC::computeInverses GetInvBQBt = " << invTimer.totalElapsedTime(),1);
279}
280
282void PresLaplaceLSCStrategy::initializeFromParameterList(const Teuchos::ParameterList & pl,const InverseLibrary & invLib)
283{
284 // get string specifying inverse
285 std::string invStr="Amesos", invVStr="", invPStr="";
286 bool useLDU = false;
287 scaleType_ = AbsRowSum;
288
289 // "parse" the parameter list
290 if(pl.isParameter("Inverse Type"))
291 invStr = pl.get<std::string>("Inverse Type");
292 if(pl.isParameter("Inverse Velocity Type"))
293 invVStr = pl.get<std::string>("Inverse Velocity Type");
294 if(pl.isParameter("Inverse Pressure Type"))
295 invPStr = pl.get<std::string>("Inverse Pressure Type");
296 if(pl.isParameter("Use LDU"))
297 useLDU = pl.get<bool>("Use LDU");
298 if(pl.isParameter("Use Mass Scaling"))
299 useMass_ = pl.get<bool>("Use Mass Scaling");
300 if(pl.isParameter("Eigen Solver Iterations"))
301 eigSolveParam_ = pl.get<int>("Eigen Solver Iterations");
302 if(pl.isParameter("Scaling Type")) {
303 scaleType_ = getDiagonalType(pl.get<std::string>("Scaling Type"));
304 TEUCHOS_TEST_FOR_EXCEPT(scaleType_==NotDiag);
305 }
306
307 // set defaults as needed
308 if(invVStr=="") invVStr = invStr;
309 if(invPStr=="") invPStr = invStr;
310
311 Teko_DEBUG_MSG_BEGIN(5)
312 DEBUG_STREAM << "LSC Inverse Strategy Parameters: " << std::endl;
313 DEBUG_STREAM << " inv v type = \"" << invVStr << "\"" << std::endl;
314 DEBUG_STREAM << " inv p type = \"" << invPStr << "\"" << std::endl;
315 DEBUG_STREAM << " use ldu = " << useLDU << std::endl;
316 DEBUG_STREAM << " use mass = " << useMass_ << std::endl;
317 DEBUG_STREAM << " scale type = " << getDiagonalName(scaleType_) << std::endl;
318 DEBUG_STREAM << "LSC Pressure Laplace Strategy Parameter list: " << std::endl;
319 pl.print(DEBUG_STREAM);
320 Teko_DEBUG_MSG_END()
321
322 // build velocity inverse factory
323 invFactoryV_ = invLib.getInverseFactory(invVStr);
324 invFactoryP_ = invFactoryV_; // by default these are the same
325 if(invVStr!=invPStr) // if different, build pressure inverse factory
326 invFactoryP_ = invLib.getInverseFactory(invPStr);
327
328 // set other parameters
329 setUseFullLDU(useLDU);
330}
331
333Teuchos::RCP<Teuchos::ParameterList> PresLaplaceLSCStrategy::getRequestedParameters() const
334{
335 Teko_DEBUG_SCOPE("PresLaplaceLSCStrategy::getRequestedParameters",10);
336 Teuchos::RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList());
337
338 // grab parameters from F solver
339 RCP<Teuchos::ParameterList> fList = invFactoryV_->getRequestedParameters();
340 if(fList!=Teuchos::null) {
341 Teuchos::ParameterList::ConstIterator itr;
342 for(itr=fList->begin();itr!=fList->end();++itr)
343 pl->setEntry(itr->first,itr->second);
344 }
345
346 // grab parameters from S solver
347 RCP<Teuchos::ParameterList> sList = invFactoryP_->getRequestedParameters();
348 if(sList!=Teuchos::null) {
349 Teuchos::ParameterList::ConstIterator itr;
350 for(itr=sList->begin();itr!=sList->end();++itr)
351 pl->setEntry(itr->first,itr->second);
352 }
353
354 // use the mass matrix
355 if(useMass_)
356 pl->set(getVelocityMassString(), false,"Velocity mass matrix");
357 pl->set(getPressureLaplaceString(), false,"Pressure Laplacian matrix");
358
359 return pl;
360}
361
363bool PresLaplaceLSCStrategy::updateRequestedParameters(const Teuchos::ParameterList & pl)
364{
365 Teko_DEBUG_SCOPE("PresLaplaceLSCStrategy::updateRequestedParameters",10);
366 bool result = true;
367
368 // update requested parameters in solvers
369 result &= invFactoryV_->updateRequestedParameters(pl);
370 result &= invFactoryP_->updateRequestedParameters(pl);
371
372 Teuchos::ParameterList hackList(pl);
373
374 // get required operator acknowledgment...user must set these to true
375 bool plo = hackList.get<bool>(getPressureLaplaceString(),false);
376
377 bool vmo = true;
378 if(useMass_)
379 vmo = hackList.get<bool>(getVelocityMassString(),false);
380
381 if(not plo) { Teko_DEBUG_MSG("User must acknowledge the use of the \"" << getPressureLaplaceString() << "\"!",0); }
382 if(not vmo) { Teko_DEBUG_MSG("User must acknowledge the use of the \"" << getVelocityMassString() << "\"!",0); }
383
384 result &= (plo & vmo);
385
386 return result;
387}
388
389} // end namespace NS
390} // 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.
@ AbsRowSum
Specifies that the diagonal entry is .
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.