Teko Version of the Day
Loading...
Searching...
No Matches
Teko_PCDStrategy.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_PCDStrategy.hpp"
48
49#include "Teuchos_TimeMonitor.hpp"
50#include "Teko_Utilities.hpp"
51
52namespace Teko {
53namespace NS {
54
55using Teuchos::TimeMonitor;
56
57Teuchos::RCP<Teuchos::Time> PCDStrategy::initTimer_;
58Teuchos::RCP<Teuchos::Time> PCDStrategy::invSTimer_;
59Teuchos::RCP<Teuchos::Time> PCDStrategy::invFTimer_;
60Teuchos::RCP<Teuchos::Time> PCDStrategy::opsTimer_;
61
63{
64 if(initTimer_==Teuchos::null)
65 initTimer_ = TimeMonitor::getNewTimer("PCDStrategy::initializePrec");
66
67 if(invSTimer_==Teuchos::null)
68 invSTimer_ = TimeMonitor::getNewTimer("PCDStrategy::initializePrec invS");
69
70 if(invFTimer_==Teuchos::null)
71 invFTimer_ = TimeMonitor::getNewTimer("PCDStrategy::initializePrec invF");
72
73 if(opsTimer_==Teuchos::null)
74 opsTimer_ = TimeMonitor::getNewTimer("PCDStrategy::initializePrec buildOps");
75}
76
77PCDStrategy::PCDStrategy() : massInverseType_(Diagonal), schurCompOrdering_(false)
78{
79 pcdParams_ = Teuchos::rcp(new Teuchos::ParameterList);
80 lapParams_ = Teuchos::rcp(new Teuchos::ParameterList);
81
82 lapParams_->set("Name",getPressureLaplaceString());
83 pcdParams_->set("Name",getPCDString());
84
86}
87
89PCDStrategy::PCDStrategy(const Teuchos::RCP<InverseFactory> & invFA,
90 const Teuchos::RCP<InverseFactory> & invS)
91 : invFactoryF_(invFA), invFactoryS_(invS), massInverseType_(Diagonal), schurCompOrdering_(false)
92{
93 pcdParams_ = Teuchos::rcp(new Teuchos::ParameterList);
94 lapParams_ = Teuchos::rcp(new Teuchos::ParameterList);
95
96 lapParams_->set("Name",getPressureLaplaceString());
97 pcdParams_->set("Name",getPCDString());
98
100}
101
103const Teko::LinearOp
104PCDStrategy::getHatInvA00(const Teko::BlockedLinearOp & A,BlockPreconditionerState & state) const
105{
106 initializeState(A,state);
107
108 return state.getModifiableOp("invF");
109}
110
112const Teko::LinearOp
113PCDStrategy::getTildeInvA00(const Teko::BlockedLinearOp & A,BlockPreconditionerState & state) const
114{
115 initializeState(A,state);
116
117 return state.getModifiableOp("invF");
118}
119
121const Teko::LinearOp
122PCDStrategy::getInvS(const Teko::BlockedLinearOp & A,BlockPreconditionerState & state) const
123{
124 initializeState(A,state);
125
126 return state.getLinearOp("invS");
127}
128
129void PCDStrategy::initializeState(const Teko::BlockedLinearOp & A,BlockPreconditionerState & state) const
130{
131 Teko_DEBUG_SCOPE("PCDStrategy::initializeState",10);
132 TEUCHOS_ASSERT(getRequestHandler()!=Teuchos::null);
133
134 std::string pcdStr = getPCDString();
135 std::string presLapStr = getPressureLaplaceString();
136 std::string presMassStr = getPressureMassString();
137
138 // no work to be done
139 if(state.isInitialized())
140 return;
141
142 Teuchos::TimeMonitor timer(*initTimer_,true);
143
144 // extract sub blocks
145 LinearOp F = Teko::getBlock(0,0,A);
146 LinearOp Bt = Teko::getBlock(0,1,A);
147 LinearOp B = Teko::getBlock(1,0,A);
148 LinearOp C = Teko::getBlock(1,1,A);
149
150 LinearOp Qp = getRequestHandler()->request<LinearOp>(presMassStr);
151 TEUCHOS_ASSERT(Qp!=Teuchos::null);
152
153 // build the inverse Laplacian complement
155 LinearOp iQp;
156 if(massInverseType_==NotDiag) {
157 ModifiableLinearOp & invMass = state.getModifiableOp("invMass");
158 Teko_DEBUG_SCOPE("Building inv(Mass)",10);
159
160 if(invMass==Teuchos::null)
161 invMass = buildInverse(*invFactoryS_,Qp);
162 else
163 rebuildInverse(*invFactoryS_,Qp,invMass);
164
165 iQp = invMass;
166 }
167 else {
168 Teko_DEBUG_MSG("Building inverse mass of type \"" << Teko::getDiagonalName(massInverseType_) << "\"",10);
169 iQp = getInvDiagonalOp(Qp,massInverseType_);
170 }
171
172 // build the inverse Laplacian complement
174 ModifiableLinearOp & invLaplace = state.getModifiableOp("invLaplace");
175 {
176 Teuchos::TimeMonitor timerInvS(*invSTimer_,true);
177
178 // LinearOp laplace = getRequestHandler()->request<Teko::LinearOp>(presLapStr);
179 LinearOp laplace = getRequestHandler()->request<Teko::LinearOp>(RequestMesg(lapParams_));
180 TEUCHOS_ASSERT(laplace!=Teuchos::null);
181 if(invLaplace==Teuchos::null)
182 invLaplace = buildInverse(*invFactoryS_,laplace);
183 else
184 rebuildInverse(*invFactoryS_,laplace,invLaplace);
185 }
186
187 // build the inverse Schur complement
189 {
190 Teko_DEBUG_SCOPE("Building S",10);
191 Teuchos::TimeMonitor timerS(*opsTimer_,true);
192
193 // build Schur-complement
194 // LinearOp pcd = getRequestHandler()->request<Teko::LinearOp>(pcdStr);
195 LinearOp pcd = getRequestHandler()->request<Teko::LinearOp>(RequestMesg(pcdParams_));
196 TEUCHOS_ASSERT(pcd!=Teuchos::null);
197 LinearOp invL = invLaplace;
198
199 LinearOp invS;
200 if(schurCompOrdering_==false)
201 invS = multiply(iQp,pcd,invL);
202 else
203 invS = multiply(invL,pcd,iQp);
204
205 state.addLinearOp("invS",invS);
206 }
207
208 // build inverse F
210 {
211 Teko_DEBUG_SCOPE("Building inv(F)",10);
212 Teuchos::TimeMonitor timerInvF(*invFTimer_,true);
213
214 ModifiableLinearOp & invF = state.getModifiableOp("invF");
215 if(invF==Teuchos::null)
216 invF = buildInverse(*invFactoryF_,F);
217 else
218 rebuildInverse(*invFactoryF_,F,invF);
219 }
220
221 // mark state as initialized
222 state.setInitialized(true);
223}
224
236void PCDStrategy::initializeFromParameterList(const Teuchos::ParameterList & pl,
237 const InverseLibrary & invLib)
238{
239 Teko_DEBUG_SCOPE("PCDStrategy::initializeFromParameterList",10);
240
241 std::string invStr="Amesos", invFStr="", invSStr="";
242 massInverseType_ = Diagonal;
243
244 // "parse" the parameter list
245 if(pl.isParameter("Inverse Type"))
246 invStr = pl.get<std::string>("Inverse Type");
247 if(pl.isParameter("Inverse F Type"))
248 invFStr = pl.get<std::string>("Inverse F Type");
249 if(pl.isParameter("Inverse Laplace Type"))
250 invSStr = pl.get<std::string>("Inverse Laplace Type");
251 if(pl.isParameter("Inverse Mass Type")) {
252 std::string massInverseStr = pl.get<std::string>("Inverse Mass Type");
253
254 // build inverse types
255 massInverseType_ = getDiagonalType(massInverseStr);
256 }
257 if(pl.isParameter("Flip Schur Complement Ordering"))
258 schurCompOrdering_ = pl.get<bool>("Flip Schur Complement Ordering");
259
260 // set defaults as needed
261 if(invFStr=="") invFStr = invStr;
262 if(invSStr=="") invSStr = invStr;
263
264 // read pressure laplace parameters
265 if(pl.isSublist("Pressure Laplace Parameters"))
266 lapParams_ = Teuchos::rcp(new Teuchos::ParameterList(pl.sublist("Pressure Laplace Parameters")));
267 else
268 lapParams_ = Teuchos::rcp(new Teuchos::ParameterList);
269
270 // read pcd operator parameters
271 if(pl.isSublist("Pressure Convection Diffusion Parameters"))
272 pcdParams_ = Teuchos::rcp(new Teuchos::ParameterList(pl.sublist("Pressure Convection Diffusion Parameters")));
273 else
274 pcdParams_ = Teuchos::rcp(new Teuchos::ParameterList);
275
276 // The user should not have already added this parameters
277 TEUCHOS_TEST_FOR_EXCEPTION(lapParams_->isParameter("Name"),std::logic_error,
278 "Teko: Parameter \"Name\" is not allowed in the sublist \""+lapParams_->name()+"\"");
279 TEUCHOS_TEST_FOR_EXCEPTION(lapParams_->isParameter("Tag"),std::logic_error,
280 "Teko: Parameter \"Tag\" is not allowed in the sublist \""+lapParams_->name()+"\"");
281 TEUCHOS_TEST_FOR_EXCEPTION(pcdParams_->isParameter("Name"),std::logic_error,
282 "Teko: Parameter \"Name\" is not allowed in the sublist \""+pcdParams_->name()+"\"");
283 TEUCHOS_TEST_FOR_EXCEPTION(pcdParams_->isParameter("Tag"),std::logic_error,
284 "Teko: Parameter \"Tag\" is not allowed in the sublist \""+pcdParams_->name()+"\"");
285
286 Teko_DEBUG_MSG_BEGIN(5)
287 DEBUG_STREAM << "PCD Strategy Parameters: " << std::endl;
288 DEBUG_STREAM << " inv type = \"" << invStr << "\"" << std::endl;
289 DEBUG_STREAM << " inv F type = \"" << invFStr << "\"" << std::endl;
290 DEBUG_STREAM << " inv Laplace type = \"" << invSStr << "\"" << std::endl;
291 DEBUG_STREAM << " inv Mass type = \"" << Teko::getDiagonalName(massInverseType_) << "\"" << std::endl;
292 DEBUG_STREAM << "PCD Strategy Parameter list: " << std::endl;
293 pl.print(DEBUG_STREAM);
294 Teko_DEBUG_MSG_END()
295
296 // build velocity inverse factory
297 invFactoryF_ = invLib.getInverseFactory(invFStr);
298
299 if(invFStr==invSStr)
300 invFactoryS_ = invFactoryF_;
301 else
302 invFactoryS_ = invLib.getInverseFactory(invSStr);
303
304 lapParams_->set("Name",getPressureLaplaceString());
305 pcdParams_->set("Name",getPCDString());
306
307 // setup a request for required operators
308 getRequestHandler()->preRequest<Teko::LinearOp>(getPressureMassString());
309 // getRequestHandler()->preRequest<Teko::LinearOp>(getPCDString());
310 // getRequestHandler()->preRequest<Teko::LinearOp>(getPressureLaplaceString());
311 getRequestHandler()->preRequest<Teko::LinearOp>(Teko::RequestMesg(lapParams_));
312 getRequestHandler()->preRequest<Teko::LinearOp>(Teko::RequestMesg(pcdParams_));
313}
314
316Teuchos::RCP<Teuchos::ParameterList> PCDStrategy::getRequestedParameters() const
317{
318 TEUCHOS_ASSERT(false);
319
320 return Teuchos::null;
321}
322
324bool PCDStrategy::updateRequestedParameters(const Teuchos::ParameterList & /* pl */)
325{
326 TEUCHOS_ASSERT(false);
327
328 return true;
329}
330
331} // end namespace NS
332} // 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.
Teuchos::RCP< RequestHandler > getRequestHandler() const
This method gets the request handler uses by this object.
virtual const Teko::LinearOp getTildeInvA00(const Teko::BlockedLinearOp &A, BlockPreconditionerState &state) const
virtual void initializeFromParameterList(const Teuchos::ParameterList &settings, const InverseLibrary &invLib)
This function builds the internals of the state from a parameter list.
virtual const Teko::LinearOp getInvS(const Teko::BlockedLinearOp &A, BlockPreconditionerState &state) const
virtual const Teko::LinearOp getHatInvA00(const Teko::BlockedLinearOp &A, BlockPreconditionerState &state) const
virtual Teuchos::RCP< Teuchos::ParameterList > getRequestedParameters() const
Request the additional parameters this preconditioner factory needs.
virtual bool updateRequestedParameters(const Teuchos::ParameterList &pl)
Update this object with the fields from a parameter list.
void initializeState(const Teko::BlockedLinearOp &A, BlockPreconditionerState &state) const
PCDStrategy()
default Constructor