47 #include "Teko_SIMPLEPreconditionerFactory.hpp"
50 #include "Teko_InverseFactory.hpp"
51 #include "Teko_BlockLowerTriInverseOp.hpp"
52 #include "Teko_BlockUpperTriInverseOp.hpp"
54 #ifdef TEKO_HAVE_EPETRA
55 #include "Teko_DiagonalPreconditionerFactory.hpp"
58 #include "Teuchos_Time.hpp"
68 : invVelFactory_(inverse),
69 invPrsFactory_(inverse),
71 fInverseType_(Diagonal),
74 SIMPLEPreconditionerFactory ::SIMPLEPreconditionerFactory(
const RCP<InverseFactory>& invVFact,
75 const RCP<InverseFactory>& invPFact,
77 : invVelFactory_(invVFact),
78 invPrsFactory_(invPFact),
80 fInverseType_(Diagonal),
83 SIMPLEPreconditionerFactory::SIMPLEPreconditionerFactory()
84 : alpha_(1.0), fInverseType_(Diagonal), useMass_(false) {}
89 Teko_DEBUG_SCOPE(
"SIMPLEPreconditionerFactory::buildPreconditionerOperator", 10);
90 Teko_DEBUG_EXPR(Teuchos::Time timer(
""));
92 int rows = blockRowCount(blockOp);
93 int cols = blockColCount(blockOp);
95 TEUCHOS_ASSERT(rows == 2);
96 TEUCHOS_ASSERT(cols == 2);
98 bool buildExplicitSchurComplement =
true;
101 const LinearOp F = getBlock(0, 0, blockOp);
102 const LinearOp Bt = getBlock(0, 1, blockOp);
103 const LinearOp B = getBlock(1, 0, blockOp);
104 const LinearOp C = getBlock(1, 1, blockOp);
108 TEUCHOS_ASSERT(massMatrix_ != Teuchos::null);
113 std::string fApproxStr =
"<error>";
115 if (fInverseType_ == NotDiag) {
117 fApproxStr = customHFactory_->toString();
120 buildExplicitSchurComplement =
false;
121 }
else if (fInverseType_ == BlkDiag) {
122 #ifdef TEKO_HAVE_EPETRA
135 buildExplicitSchurComplement =
true;
138 throw std::logic_error(
139 "SIMPLEPreconditionerFactory fInverseType_ == "
140 "BlkDiag but EPETRA is turned off!");
145 H = getInvDiagonalOp(matF, fInverseType_);
146 fApproxStr = getDiagonalName(fInverseType_);
151 RCP<const Teuchos::ParameterList> pl = state.getParameterList();
153 if (pl->isParameter(
"stepsize")) {
155 double stepsize = pl->get<
double>(
"stepsize");
158 if (stepsize > 0.0) H = scale(stepsize, H);
165 if (buildExplicitSchurComplement) {
171 mHBt = explicitMultiply(H, Bt, mHBt);
175 BHBt = explicitMultiply(B, HBt, BHBt);
178 mhatS = explicitAdd(C, scale(-1.0, BHBt), mhatS);
182 HBt = multiply(H, Bt);
184 hatS = add(C, scale(-1.0, multiply(B, HBt)));
187 Teko::ModifiableLinearOp& precInvF = state.
getModifiableOp(
"precInvF");
189 if(precInvF == Teuchos::null){
190 precInvF = precVelFactory_->buildInverse(F);
199 if (invF == Teuchos::null){
200 if(precInvF.is_null()){
207 if(precInvF.is_null()){
214 Teko::ModifiableLinearOp& precInvS = state.
getModifiableOp(
"precInvS");
216 if(precInvS == Teuchos::null){
217 precInvS = precPrsFactory_->buildInverse(hatS);
226 if (invS == Teuchos::null){
227 if(precInvS == Teuchos::null){
234 if(precInvS == Teuchos::null){
241 std::vector<LinearOp> invDiag(2);
244 BlockedLinearOp L = zeroBlockedOp(blockOp);
245 setBlock(1, 0, L, B);
250 LinearOp invL = createBlockLowerTriInverseOp(L, invDiag);
253 BlockedLinearOp U = zeroBlockedOp(blockOp);
254 setBlock(0, 1, U, scale(1.0 / alpha_, HBt));
257 invDiag[0] = identity(rangeSpace(invF));
258 invDiag[1] = scale(alpha_, identity(rangeSpace(invS)));
259 LinearOp invU = createBlockUpperTriInverseOp(U, invDiag);
262 return multiply(invU, invL,
"SIMPLE_" + fApproxStr);
271 customHFactory_ = Teuchos::null;
272 fInverseType_ = Diagonal;
275 std::string invStr =
"", invVStr =
"", invPStr =
"", precVStr =
"", precPStr =
"";
279 if (pl.isParameter(
"Inverse Type")) invStr = pl.get<std::string>(
"Inverse Type");
280 if (pl.isParameter(
"Inverse Velocity Type"))
281 invVStr = pl.get<std::string>(
"Inverse Velocity Type");
282 if (pl.isParameter(
"Preconditioner Velocity Type"))
283 precVStr = pl.get<std::string>(
"Preconditioner Velocity Type");
284 if (pl.isParameter(
"Inverse Pressure Type"))
285 invPStr = pl.get<std::string>(
"Inverse Pressure Type");
286 if (pl.isParameter(
"Preconditioner Pressure Type"))
287 precPStr = pl.get<std::string>(
"Preconditioner Pressure Type");
288 if (pl.isParameter(
"Alpha")) alpha_ = pl.get<
double>(
"Alpha");
289 if (pl.isParameter(
"Explicit Velocity Inverse Type")) {
290 std::string fInverseStr = pl.get<std::string>(
"Explicit Velocity Inverse Type");
293 fInverseType_ = getDiagonalType(fInverseStr);
294 if (fInverseType_ == NotDiag) customHFactory_ = invLib->getInverseFactory(fInverseStr);
297 if (fInverseType_ == BlkDiag) BlkDiagList_ = pl.sublist(
"H options");
299 if (pl.isParameter(
"Use Mass Scaling")) useMass_ = pl.get<
bool>(
"Use Mass Scaling");
301 Teko_DEBUG_MSG_BEGIN(5) DEBUG_STREAM <<
"SIMPLE Parameters: " << std::endl;
302 DEBUG_STREAM <<
" inv type = \"" << invStr <<
"\"" << std::endl;
303 DEBUG_STREAM <<
" inv v type = \"" << invVStr <<
"\"" << std::endl;
304 DEBUG_STREAM <<
" prec v type = \"" << precVStr <<
"\"" << std::endl;
305 DEBUG_STREAM <<
" inv p type = \"" << invPStr <<
"\"" << std::endl;
306 DEBUG_STREAM <<
" prec p type = \"" << precPStr <<
"\"" << std::endl;
307 DEBUG_STREAM <<
" alpha = " << alpha_ << std::endl;
308 DEBUG_STREAM <<
" use mass = " << useMass_ << std::endl;
309 DEBUG_STREAM <<
" vel scaling = " << getDiagonalName(fInverseType_) << std::endl;
310 DEBUG_STREAM <<
"SIMPLE Parameter list: " << std::endl;
311 pl.print(DEBUG_STREAM);
315 #if defined(Teko_ENABLE_Amesos)
316 if (invStr ==
"") invStr =
"Amesos";
317 #elif defined(Teko_ENABLE_Amesos2)
318 if (invStr ==
"") invStr =
"Amesos2";
321 if (invVStr ==
"") invVStr = invStr;
322 if (invPStr ==
"") invPStr = invStr;
325 RCP<InverseFactory> invVFact, invPFact;
328 invVFact = invLib->getInverseFactory(invVStr);
330 if (invVStr != invPStr)
331 invPFact = invLib->getInverseFactory(invPStr);
333 RCP<InverseFactory> precVFact, precPFact;
335 precVFact = invLib->getInverseFactory(precVStr);
338 precPFact = invLib->getInverseFactory(precPStr);
341 invVelFactory_ = invVFact;
342 invPrsFactory_ = invPFact;
344 precVelFactory_ = precVFact;
345 precPrsFactory_ = precPFact;
349 rh->preRequest<Teko::LinearOp>(Teko::RequestMesg(
"Velocity Mass Matrix"));
350 Teko::LinearOp mass = rh->request<Teko::LinearOp>(Teko::RequestMesg(
"Velocity Mass Matrix"));
357 Teuchos::RCP<Teuchos::ParameterList> result;
358 Teuchos::RCP<Teuchos::ParameterList> pl = rcp(
new Teuchos::ParameterList());
362 RCP<Teuchos::ParameterList> vList = invVelFactory_->getRequestedParameters();
363 if (vList != Teuchos::null) {
364 Teuchos::ParameterList::ConstIterator itr;
365 for (itr = vList->begin(); itr != vList->end(); ++itr) pl->setEntry(itr->first, itr->second);
370 if (precVelFactory_ != Teuchos::null) {
371 RCP<Teuchos::ParameterList> vList = precVelFactory_->getRequestedParameters();
372 if (vList != Teuchos::null) {
373 Teuchos::ParameterList::ConstIterator itr;
374 for (itr = vList->begin(); itr != vList->end(); ++itr) pl->setEntry(itr->first, itr->second);
381 RCP<Teuchos::ParameterList> pList = invPrsFactory_->getRequestedParameters();
382 if (pList != Teuchos::null) {
383 Teuchos::ParameterList::ConstIterator itr;
384 for (itr = pList->begin(); itr != pList->end(); ++itr) pl->setEntry(itr->first, itr->second);
389 if (precPrsFactory_ != Teuchos::null) {
390 RCP<Teuchos::ParameterList> pList = precPrsFactory_->getRequestedParameters();
391 if (pList != Teuchos::null) {
392 Teuchos::ParameterList::ConstIterator itr;
393 for (itr = pList->begin(); itr != pList->end(); ++itr) pl->setEntry(itr->first, itr->second);
399 if (customHFactory_ != Teuchos::null) {
400 RCP<Teuchos::ParameterList> hList = customHFactory_->getRequestedParameters();
401 if (hList != Teuchos::null) {
402 Teuchos::ParameterList::ConstIterator itr;
403 for (itr = hList->begin(); itr != hList->end(); ++itr) pl->setEntry(itr->first, itr->second);
413 Teko_DEBUG_SCOPE(
"InvLSCStrategy::updateRequestedParameters", 10);
417 result &= invVelFactory_->updateRequestedParameters(pl);
418 result &= invPrsFactory_->updateRequestedParameters(pl);
420 result &= precVelFactory_->updateRequestedParameters(pl);
422 result &= precPrsFactory_->updateRequestedParameters(pl);
423 if (customHFactory_ != Teuchos::null) result &= customHFactory_->updateRequestedParameters(pl);
void rebuildInverse(const InverseFactory &factory, const LinearOp &A, InverseLinearOp &invA)
virtual bool updateRequestedParameters(const Teuchos::ParameterList &pl)
For assisting in construction of the preconditioner.
Preconditioner factory for building explcit inverse of diagonal operators. This includes block operat...
An implementation of a state object for block preconditioners.
virtual void setMassMatrix(Teko::LinearOp &mass)
Set the mass matrix for this factory.
virtual void addModifiableOp(const std::string &name, const Teko::ModifiableLinearOp &mlo)
Add a named operator to the state object.
InverseLinearOp buildInverse(const InverseFactory &factory, const LinearOp &A)
Build an inverse operator using a factory and a linear operator.
virtual void initializeFromParameterList(const Teuchos::ParameterList &pl)
Initialize from a parameter list.
virtual void initializeFromParameterList(const Teuchos::ParameterList &pl)
Initialize from a parameter list.
Teuchos::RCP< const InverseLibrary > getInverseLibrary() const
Get the inverse library used by this preconditioner factory.
virtual Teuchos::RCP< Teuchos::ParameterList > getRequestedParameters() const
For assisting in construction of the preconditioner.
LinearOp buildPreconditionerOperator(BlockedLinearOp &blo, BlockPreconditionerState &state) const
Function that is called to build the preconditioner for the linear operator that is passed in...
Teuchos::RCP< RequestHandler > getRequestHandler() const
Get the request handler with pointers to the appropriate callbacks.
SIMPLEPreconditionerFactory()
Default constructor.
LinearOp buildPreconditionerOperator(LinearOp &lo, PreconditionerState &state) const
virtual Teko::ModifiableLinearOp & getModifiableOp(const std::string &name)
Add a named operator to the state object.