102LinearOp TimingsSIMPLEPreconditionerFactory
103 ::buildPreconditionerOperator(BlockedLinearOp & blockOp,
106 constrTotal_.start();
108 Teko_DEBUG_SCOPE(
"TimingsSIMPLEPreconditionerFactory::buildPreconditionerOperator",10);
109 Teko_DEBUG_EXPR(Teuchos::Time timer(
""));
114 TEUCHOS_ASSERT(rows==2);
115 TEUCHOS_ASSERT(cols==2);
117 bool buildExplicitSchurComplement =
true;
120 const LinearOp F =
getBlock(0,0,blockOp);
121 const LinearOp Bt =
getBlock(0,1,blockOp);
122 const LinearOp B =
getBlock(1,0,blockOp);
123 const LinearOp C =
getBlock(1,1,blockOp);
127 TEUCHOS_ASSERT(massMatrix_!=Teuchos::null);
132 std::string fApproxStr =
"<error>";
136 fApproxStr = customHFactory_->toString();
139 buildExplicitSchurComplement =
false;
141 else if(fInverseType_==
BlkDiag) {
142#ifdef TEKO_HAVE_EPETRA
146 Hfact.initializeFromParameterList(BlkDiagList_);
147 H = Hfact.buildPreconditionerOperator(matF,Hstate);
149 buildExplicitSchurComplement =
true;
152 throw std::logic_error(
"TimingsSIMPLEPreconditionerFactory fInverseType_ == "
153 "BlkDiag but EPETRA is turned off!");
159 H = getInvDiagonalOp(matF,fInverseType_);
161 fApproxStr = getDiagonalName(fInverseType_);
166 RCP<const Teuchos::ParameterList> pl = state.getParameterList();
168 if(pl->isParameter(
"stepsize")) {
170 double stepsize = pl->get<
double>(
"stepsize");
174 H = scale(stepsize,H);
181 if(buildExplicitSchurComplement) {
182 ModifiableLinearOp & mHBt = state.getModifiableOp(
"HBt");
183 ModifiableLinearOp & mhatS = state.getModifiableOp(
"hatS");
184 ModifiableLinearOp & BHBt = state.getModifiableOp(
"BHBt");
188 mHBt = explicitMultiply(H,Bt,mHBt);
194 BHBt = explicitMultiply(B,HBt,BHBt);
199 mhatS = explicitAdd(C,scale(-1.0,BHBt),mhatS);
205 HBt = multiply(H,Bt);
207 hatS = add(C,scale(-1.0,multiply(B,HBt)));
211 if(timed_HBt_==Teuchos::null) {
215 timed_HBt_->setLinearOp(HBt);
219 if(timed_B_==Teuchos::null) {
223 timed_B_->setLinearOp(B);
227 ModifiableLinearOp & invF = state.getModifiableOp(
"invF");
229 if(invF==Teuchos::null) {
237 timed_invF_->setLinearOp(invF);
242 ModifiableLinearOp & invS = state.getModifiableOp(
"invS");
244 if(invS==Teuchos::null) {
252 timed_invS_->setLinearOp(invS);
256 std::vector<LinearOp> invDiag(2);
259 BlockedLinearOp L = zeroBlockedOp(blockOp);
263 invDiag[0] = timed_invF_;
264 invDiag[1] = timed_invS_;
265 LinearOp invL = createBlockLowerTriInverseOp(L,invDiag);
268 BlockedLinearOp U = zeroBlockedOp(blockOp);
269 setBlock(0,1,U,scale<double>(1.0/alpha_,timed_HBt_.getConst()));
273 invDiag[1] = scale(alpha_,identity(
rangeSpace(invS)));
274 LinearOp invU = createBlockUpperTriInverseOp(U,invDiag);
277 Teko::LinearOp iU_t_iL = multiply(invU,invL,
"SIMPLE_"+fApproxStr);
280 if(timed_iU_t_iL_==Teuchos::null)
281 timed_iU_t_iL_ = Teuchos::rcp(
new DiagnosticLinearOp(getOutputStream(),iU_t_iL,
"iU_t_iL"));
283 timed_iU_t_iL_->setLinearOp(iU_t_iL);
289 return timed_iU_t_iL_;
299 customHFactory_ = Teuchos::null;
303 std::string invStr=
"", invVStr=
"", invPStr=
"";
307 if(pl.isParameter(
"Inverse Type"))
308 invStr = pl.get<std::string>(
"Inverse Type");
309 if(pl.isParameter(
"Inverse Velocity Type"))
310 invVStr = pl.get<std::string>(
"Inverse Velocity Type");
311 if(pl.isParameter(
"Inverse Pressure Type"))
312 invPStr = pl.get<std::string>(
"Inverse Pressure Type");
313 if(pl.isParameter(
"Alpha"))
314 alpha_ = pl.get<
double>(
"Alpha");
315 if(pl.isParameter(
"Explicit Velocity Inverse Type")) {
316 std::string fInverseStr = pl.get<std::string>(
"Explicit Velocity Inverse Type");
319 fInverseType_ = getDiagonalType(fInverseStr);
321 customHFactory_ = invLib->getInverseFactory(fInverseStr);
325 BlkDiagList_=pl.sublist(
"H options");
327 if(pl.isParameter(
"Use Mass Scaling"))
328 useMass_ = pl.get<
bool>(
"Use Mass Scaling");
330 Teko_DEBUG_MSG_BEGIN(5)
331 DEBUG_STREAM <<
"SIMPLE Parameters: " << std::endl;
332 DEBUG_STREAM <<
" inv type = \"" << invStr <<
"\"" << std::endl;
333 DEBUG_STREAM <<
" inv v type = \"" << invVStr <<
"\"" << std::endl;
334 DEBUG_STREAM <<
" inv p type = \"" << invPStr <<
"\"" << std::endl;
335 DEBUG_STREAM <<
" alpha = " << alpha_ << std::endl;
336 DEBUG_STREAM <<
" use mass = " << useMass_ << std::endl;
337 DEBUG_STREAM <<
" vel scaling = " << getDiagonalName(fInverseType_) << std::endl;
338 DEBUG_STREAM <<
"SIMPLE Parameter list: " << std::endl;
339 pl.print(DEBUG_STREAM);
343 if(invStr==
"") invStr =
"Amesos";
344 if(invVStr==
"") invVStr = invStr;
345 if(invPStr==
"") invPStr = invStr;
348 RCP<InverseFactory> invVFact, invPFact;
351 invVFact = invLib->getInverseFactory(invVStr);
354 invPFact = invLib->getInverseFactory(invPStr);
357 invVelFactory_ = invVFact;
358 invPrsFactory_ = invPFact;
362 rh->preRequest<Teko::LinearOp>(Teko::RequestMesg(
"Velocity Mass Matrix"));
364 = rh->request<Teko::LinearOp>(Teko::RequestMesg(
"Velocity Mass Matrix"));
372 Teuchos::RCP<Teuchos::ParameterList> result;
373 Teuchos::RCP<Teuchos::ParameterList> pl = rcp(
new Teuchos::ParameterList());
376 RCP<Teuchos::ParameterList> vList = invVelFactory_->getRequestedParameters();
377 if(vList!=Teuchos::null) {
378 Teuchos::ParameterList::ConstIterator itr;
379 for(itr=vList->begin();itr!=vList->end();++itr)
380 pl->setEntry(itr->first,itr->second);
385 RCP<Teuchos::ParameterList> pList = invPrsFactory_->getRequestedParameters();
386 if(pList!=Teuchos::null) {
387 Teuchos::ParameterList::ConstIterator itr;
388 for(itr=pList->begin();itr!=pList->end();++itr)
389 pl->setEntry(itr->first,itr->second);
394 if(customHFactory_!=Teuchos::null) {
395 RCP<Teuchos::ParameterList> hList = customHFactory_->getRequestedParameters();
396 if(hList!=Teuchos::null) {
397 Teuchos::ParameterList::ConstIterator itr;
398 for(itr=hList->begin();itr!=hList->end();++itr)
399 pl->setEntry(itr->first,itr->second);