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)));
189 if (invF == Teuchos::null)
196 if (invS == Teuchos::null)
201 std::vector<LinearOp> invDiag(2);
204 BlockedLinearOp L = zeroBlockedOp(blockOp);
205 setBlock(1, 0, L, B);
210 LinearOp invL = createBlockLowerTriInverseOp(L, invDiag);
213 BlockedLinearOp U = zeroBlockedOp(blockOp);
214 setBlock(0, 1, U, scale(1.0 / alpha_, HBt));
217 invDiag[0] = identity(rangeSpace(invF));
218 invDiag[1] = scale(alpha_, identity(rangeSpace(invS)));
219 LinearOp invU = createBlockUpperTriInverseOp(U, invDiag);
222 return multiply(invU, invL,
"SIMPLE_" + fApproxStr);
231 customHFactory_ = Teuchos::null;
232 fInverseType_ = Diagonal;
235 std::string invStr =
"", invVStr =
"", invPStr =
"";
239 if (pl.isParameter(
"Inverse Type")) invStr = pl.get<std::string>(
"Inverse Type");
240 if (pl.isParameter(
"Inverse Velocity Type"))
241 invVStr = pl.get<std::string>(
"Inverse Velocity Type");
242 if (pl.isParameter(
"Inverse Pressure Type"))
243 invPStr = pl.get<std::string>(
"Inverse Pressure Type");
244 if (pl.isParameter(
"Alpha")) alpha_ = pl.get<
double>(
"Alpha");
245 if (pl.isParameter(
"Explicit Velocity Inverse Type")) {
246 std::string fInverseStr = pl.get<std::string>(
"Explicit Velocity Inverse Type");
249 fInverseType_ = getDiagonalType(fInverseStr);
250 if (fInverseType_ == NotDiag) customHFactory_ = invLib->getInverseFactory(fInverseStr);
253 if (fInverseType_ == BlkDiag) BlkDiagList_ = pl.sublist(
"H options");
255 if (pl.isParameter(
"Use Mass Scaling")) useMass_ = pl.get<
bool>(
"Use Mass Scaling");
257 Teko_DEBUG_MSG_BEGIN(5) DEBUG_STREAM <<
"SIMPLE Parameters: " << std::endl;
258 DEBUG_STREAM <<
" inv type = \"" << invStr <<
"\"" << std::endl;
259 DEBUG_STREAM <<
" inv v type = \"" << invVStr <<
"\"" << std::endl;
260 DEBUG_STREAM <<
" inv p type = \"" << invPStr <<
"\"" << std::endl;
261 DEBUG_STREAM <<
" alpha = " << alpha_ << std::endl;
262 DEBUG_STREAM <<
" use mass = " << useMass_ << std::endl;
263 DEBUG_STREAM <<
" vel scaling = " << getDiagonalName(fInverseType_) << std::endl;
264 DEBUG_STREAM <<
"SIMPLE Parameter list: " << std::endl;
265 pl.print(DEBUG_STREAM);
269 if (invStr ==
"") invStr =
"Amesos";
270 if (invVStr ==
"") invVStr = invStr;
271 if (invPStr ==
"") invPStr = invStr;
274 RCP<InverseFactory> invVFact, invPFact;
277 invVFact = invLib->getInverseFactory(invVStr);
279 if (invVStr != invPStr)
280 invPFact = invLib->getInverseFactory(invPStr);
283 invVelFactory_ = invVFact;
284 invPrsFactory_ = invPFact;
288 rh->preRequest<Teko::LinearOp>(Teko::RequestMesg(
"Velocity Mass Matrix"));
289 Teko::LinearOp mass = rh->request<Teko::LinearOp>(Teko::RequestMesg(
"Velocity Mass Matrix"));
296 Teuchos::RCP<Teuchos::ParameterList> result;
297 Teuchos::RCP<Teuchos::ParameterList> pl = rcp(
new Teuchos::ParameterList());
300 RCP<Teuchos::ParameterList> vList = invVelFactory_->getRequestedParameters();
301 if (vList != Teuchos::null) {
302 Teuchos::ParameterList::ConstIterator itr;
303 for (itr = vList->begin(); itr != vList->end(); ++itr) pl->setEntry(itr->first, itr->second);
308 RCP<Teuchos::ParameterList> pList = invPrsFactory_->getRequestedParameters();
309 if (pList != Teuchos::null) {
310 Teuchos::ParameterList::ConstIterator itr;
311 for (itr = pList->begin(); itr != pList->end(); ++itr) pl->setEntry(itr->first, itr->second);
316 if (customHFactory_ != Teuchos::null) {
317 RCP<Teuchos::ParameterList> hList = customHFactory_->getRequestedParameters();
318 if (hList != Teuchos::null) {
319 Teuchos::ParameterList::ConstIterator itr;
320 for (itr = hList->begin(); itr != hList->end(); ++itr) pl->setEntry(itr->first, itr->second);
330 Teko_DEBUG_SCOPE(
"InvLSCStrategy::updateRequestedParameters", 10);
334 result &= invVelFactory_->updateRequestedParameters(pl);
335 result &= invPrsFactory_->updateRequestedParameters(pl);
336 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.
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.