47#include "NS/Teko_PresLaplaceLSCStrategy.hpp"
49#include "Thyra_DefaultDiagonalLinearOp.hpp"
51#include "Teuchos_Time.hpp"
52#include "Teuchos_TimeMonitor.hpp"
56#include "NS/Teko_LSCPreconditionerFactory.hpp"
59using Teuchos::rcp_dynamic_cast;
60using Teuchos::rcp_const_cast;
71PresLaplaceLSCStrategy::PresLaplaceLSCStrategy()
72 : invFactoryV_(Teuchos::null), invFactoryP_(Teuchos::null)
73 , eigSolveParam_(5), useFullLDU_(false), useMass_(false), scaleType_(
AbsRowSum)
76PresLaplaceLSCStrategy::PresLaplaceLSCStrategy(
const Teuchos::RCP<InverseFactory> & factory)
77 : invFactoryV_(factory), invFactoryP_(factory)
78 , eigSolveParam_(5), useFullLDU_(false), useMass_(false), scaleType_(
AbsRowSum)
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)
91 Teko_DEBUG_SCOPE(
"PresLaplaceLSCStrategy::buildState",10);
94 TEUCHOS_ASSERT(lscState!=0);
97 if(not lscState->isInitialized()) {
98 Teko_DEBUG_EXPR(Teuchos::Time timer(
""));
102 Teko_DEBUG_SCOPE(
"PL-LSC::buildState constructing operators",1);
103 Teko_DEBUG_EXPR(timer.start(
true));
105 initializeState(A,lscState);
107 Teko_DEBUG_EXPR(timer.stop());
108 Teko_DEBUG_MSG(
"PL-LSC::buildState BuildOpsTime = " << timer.totalElapsedTime(),1);
113 Teko_DEBUG_SCOPE(
"PL-LSC::buildState calculating inverses",1);
114 Teko_DEBUG_EXPR(timer.start(
true));
116 computeInverses(A,lscState);
118 Teko_DEBUG_EXPR(timer.stop());
119 Teko_DEBUG_MSG(
"PL-LSC::buildState BuildInvTime = " << timer.totalElapsedTime(),1);
127 return state.getModifiableOp(
"invPresLap");
132 return state.getModifiableOp(
"invPresLap");
137 return state.getModifiableOp(
"invF");
144 TEUCHOS_ASSERT(lscState!=0);
145 TEUCHOS_ASSERT(lscState->isInitialized())
147 return lscState->aiD_;
153 TEUCHOS_ASSERT(lscState!=0);
154 TEUCHOS_ASSERT(lscState->isInitialized())
156 return lscState->invMass_;
161 return getInvMass(A,state);
165void PresLaplaceLSCStrategy::initializeState(
const BlockedLinearOp & A,
LSCPrecondState * state)
const
167 Teko_DEBUG_SCOPE(
"PresLaplaceLSCStrategy::initializeState",10);
169 std::string velMassStr = getVelocityMassString();
172 const LinearOp Bt =
getBlock(0,1,A);
179 bool isStabilized = (not isZeroOp(C));
182 LinearOp massMatrix = state->getLinearOp(velMassStr);
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_);
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_);
203 if(not isStabilized) {
204 state->aiD_ = Teuchos::null;
206 state->setInitialized(
true);
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;
218 MultiVector vec_D = getDiagonal(B_idF_Bt);
219 update(-1.0,getDiagonal(C),1.0,vec_D);
220 const LinearOp invD = buildInvDiagonal(vec_D,
"inv(D)");
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);
229 Teko_DEBUG_MSG_BEGIN(5)
230 DEBUG_STREAM <<
"PL-LSC Alpha Parameter = " << state->alpha_ << std::endl;
233 state->setInitialized(
true);
241void PresLaplaceLSCStrategy::computeInverses(
const BlockedLinearOp & A,
LSCPrecondState * state)
const
243 Teko_DEBUG_SCOPE(
"PresLaplaceLSCStrategy::computeInverses",10);
244 Teko_DEBUG_EXPR(Teuchos::Time invTimer(
""));
246 std::string presLapStr = getPressureLaplaceString();
249 const LinearOp presLap = state->getLinearOp(presLapStr);
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) {
262 Teko_DEBUG_EXPR(invTimer.stop());
263 Teko_DEBUG_MSG(
"PL-LSC::computeInverses GetInvF = " << invTimer.totalElapsedTime(),1);
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) {
277 Teko_DEBUG_EXPR(invTimer.stop());
278 Teko_DEBUG_MSG(
"PL-LSC::computeInverses GetInvBQBt = " << invTimer.totalElapsedTime(),1);
282void PresLaplaceLSCStrategy::initializeFromParameterList(
const Teuchos::ParameterList & pl,
const InverseLibrary & invLib)
285 std::string invStr=
"Amesos", invVStr=
"", invPStr=
"";
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);
308 if(invVStr==
"") invVStr = invStr;
309 if(invPStr==
"") invPStr = invStr;
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);
323 invFactoryV_ = invLib.getInverseFactory(invVStr);
324 invFactoryP_ = invFactoryV_;
326 invFactoryP_ = invLib.getInverseFactory(invPStr);
329 setUseFullLDU(useLDU);
333Teuchos::RCP<Teuchos::ParameterList> PresLaplaceLSCStrategy::getRequestedParameters()
const
335 Teko_DEBUG_SCOPE(
"PresLaplaceLSCStrategy::getRequestedParameters",10);
336 Teuchos::RCP<Teuchos::ParameterList> pl = rcp(
new Teuchos::ParameterList());
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);
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);
356 pl->set(getVelocityMassString(),
false,
"Velocity mass matrix");
357 pl->set(getPressureLaplaceString(),
false,
"Pressure Laplacian matrix");
363bool PresLaplaceLSCStrategy::updateRequestedParameters(
const Teuchos::ParameterList & pl)
365 Teko_DEBUG_SCOPE(
"PresLaplaceLSCStrategy::updateRequestedParameters",10);
369 result &= invFactoryV_->updateRequestedParameters(pl);
370 result &= invFactoryP_->updateRequestedParameters(pl);
372 Teuchos::ParameterList hackList(pl);
375 bool plo = hackList.get<
bool>(getPressureLaplaceString(),
false);
379 vmo = hackList.get<
bool>(getVelocityMassString(),
false);
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); }
384 result &= (plo & vmo);
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.