47#ifndef MUELU_SIMPLESMOOTHER_DEF_HPP_
48#define MUELU_SIMPLESMOOTHER_DEF_HPP_
50#include "Teuchos_ArrayViewDecl.hpp"
51#include "Teuchos_ScalarTraits.hpp"
55#include <Xpetra_Matrix.hpp>
56#include <Xpetra_CrsMatrixWrap.hpp>
57#include <Xpetra_BlockedCrsMatrix.hpp>
58#include <Xpetra_MultiVectorFactory.hpp>
59#include <Xpetra_VectorFactory.hpp>
60#include <Xpetra_ReorderedBlockedCrsMatrix.hpp>
64#include "MueLu_Utilities.hpp"
66#include "MueLu_HierarchyUtils.hpp"
68#include "MueLu_SubBlockAFactory.hpp"
71#include "MueLu_SchurComplementFactory.hpp"
72#include "MueLu_DirectSolver.hpp"
73#include "MueLu_SmootherFactory.hpp"
74#include "MueLu_FactoryManager.hpp"
78 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
80 : type_(
"SIMPLE"), A_(
Teuchos::null) {}
82 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
85 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
87 RCP<ParameterList> validParamList = rcp(
new ParameterList());
89 validParamList->set<RCP<const FactoryBase>>(
"A", Teuchos::null,
"Generating factory of the matrix A");
90 validParamList->set<
Scalar>(
"Damping factor", 1.0,
"Damping/Scaling factor in SIMPLE");
91 validParamList->set<
LocalOrdinal>(
"Sweeps", 1,
"Number of SIMPLE sweeps (default = 1)");
92 validParamList->set<
bool>(
"UseSIMPLEC",
false,
"Use SIMPLEC instead of SIMPLE (default = false)");
94 return validParamList;
98 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
100 TEUCHOS_TEST_FOR_EXCEPTION(pos < 0,
Exceptions::RuntimeError,
"MueLu::SimpleSmoother::AddFactoryManager: parameter \'pos\' must not be negative! error.");
102 size_t myPos = Teuchos::as<size_t>(pos);
104 if (myPos < FactManager_.size()) {
106 FactManager_.at(myPos) = FactManager;
107 }
else if( myPos == FactManager_.size()) {
109 FactManager_.push_back(FactManager);
111 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
112 *out <<
"Warning: cannot add new FactoryManager at proper position " << pos <<
". The FactoryManager is just appended to the end. Check this!" << std::endl;
115 FactManager_.push_back(FactManager);
121 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
123 AddFactoryManager(FactManager, 0);
126 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
128 AddFactoryManager(FactManager, 1);
131 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
133 currentLevel.
DeclareInput(
"A",this->GetFactory(
"A").get());
135 TEUCHOS_TEST_FOR_EXCEPTION(FactManager_.size() != 2,
Exceptions::RuntimeError,
"MueLu::SimpleSmoother::DeclareInput: You have to declare two FactoryManagers with a \"Smoother\" object: One for predicting the primary variable and one for the SchurComplement system. The smoother for the SchurComplement system needs a SchurComplementFactory as input for variable \"A\". make sure that you use the same proper damping factors for omega both in the SchurComplementFactory and in the SIMPLE smoother!");
138 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
139 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
143 currentLevel.
DeclareInput(
"PreSmoother",(*it)->GetFactory(
"Smoother").get());
147 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
150 FactoryMonitor m(*
this,
"Setup blocked SIMPLE Smoother", currentLevel);
153 this->GetOStream(
Warnings0) <<
"MueLu::SimpleSmoother::Setup(): Setup() has already been called";
156 A_ = Factory::Get<RCP<Matrix> > (currentLevel,
"A");
158 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
159 TEUCHOS_TEST_FOR_EXCEPTION(bA == Teuchos::null,
Exceptions::BadCast,
"MueLu::SimpleSmoother::Setup: input matrix A is not of type BlockedCrsMatrix! error.");
162 rangeMapExtractor_ = bA->getRangeMapExtractor();
163 domainMapExtractor_ = bA->getDomainMapExtractor();
166 F_ = bA->getMatrix(0, 0);
167 G_ = bA->getMatrix(0, 1);
168 D_ = bA->getMatrix(1, 0);
169 Z_ = bA->getMatrix(1, 1);
172 bool bSIMPLEC = pL.get<
bool>(
"UseSIMPLEC");
176 RCP<Vector> diagFVector = VectorFactory::Build(F_->getRowMap());
178 F_->getLocalDiagCopy(*diagFVector);
199 RCP<BlockedVector> bdiagFinv = Teuchos::rcp_dynamic_cast<BlockedVector>(diagFinv_);
200 if(bdiagFinv.is_null() ==
false && bdiagFinv->getBlockedMap()->getNumMaps() == 1) {
201 RCP<Vector> nestedVec = bdiagFinv->getMultiVector(0,bdiagFinv->getBlockedMap()->getThyraMode())->getVectorNonConst(0);
202 diagFinv_.swap(nestedVec);
208 RCP<const FactoryManagerBase> velpredictFactManager = FactManager_.at(0);
210 velPredictSmoo_ = currentLevel.
Get< RCP<SmootherBase> >(
"PreSmoother",velpredictFactManager->GetFactory(
"Smoother").get());
213 RCP<const FactoryManagerBase> schurFactManager = FactManager_.at(1);
215 schurCompSmoo_ = currentLevel.
Get< RCP<SmootherBase> >(
"PreSmoother", schurFactManager->GetFactory(
"Smoother").get());
221 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
225 "MueLu::SimpleSmoother::Apply(): Setup() has not been called");
228 RCP<MultiVector> rcpDebugX = Teuchos::rcpFromRef(X);
229 RCP<const MultiVector> rcpDebugB = Teuchos::rcpFromRef(B);
230 RCP<BlockedMultiVector> rcpBDebugX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpDebugX);
231 RCP<const BlockedMultiVector> rcpBDebugB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpDebugB);
232 if(rcpBDebugB.is_null() ==
false) {
237 if(rcpBDebugX.is_null() ==
false) {
244 const SC zero = Teuchos::ScalarTraits<SC>::zero();
245 const SC one = Teuchos::ScalarTraits<SC>::one();
253 bool bRangeThyraMode = rangeMapExtractor_->getThyraMode();
254 bool bDomainThyraMode = domainMapExtractor_->getThyraMode();
259 RCP<MultiVector> rcpX = Teuchos::rcpFromRef(X);
260 RCP<const MultiVector> rcpB = Teuchos::rcpFromRef(B);
263 bool bCopyResultX =
false;
264 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
266 "MueLu::BlockedGaussSeidelSmoother::Apply(): A_ must be a BlockedCrsMatrix");
267 RCP<BlockedMultiVector> bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
268 RCP<const BlockedMultiVector> bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
270 if(bX.is_null() ==
true) {
271 RCP<MultiVector> test = Teuchos::rcp(
new BlockedMultiVector(bA->getBlockedDomainMap(),rcpX));
276 if(bB.is_null() ==
true) {
277 RCP<const MultiVector> test = Teuchos::rcp(
new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB));
282 bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
283 bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
286 RCP<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rbA =
287 Teuchos::rcp_dynamic_cast<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(bA);
288 if(rbA.is_null() ==
false) {
290 Teuchos::RCP<const Xpetra::BlockReorderManager > brm = rbA->getBlockReorderManager();
293 if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps()) {
295 Teuchos::RCP<MultiVector> test = buildReorderedBlockedMultiVector(brm, bX);
298 if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps()) {
300 Teuchos::RCP<const MultiVector> test = buildReorderedBlockedMultiVector(brm, bB);
309 RCP<MultiVector> residual = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors());
310 RCP<BlockedMultiVector> bresidual = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(residual);
311 RCP<MultiVector> r1 = bresidual->getMultiVector(0,bRangeThyraMode);
312 RCP<MultiVector> r2 = bresidual->getMultiVector(1,bRangeThyraMode);
315 RCP<MultiVector> xtilde = MultiVectorFactory::Build(rcpX->getMap(), rcpX->getNumVectors());
316 RCP<BlockedMultiVector> bxtilde = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(xtilde);
317 RCP<MultiVector> xtilde1 = bxtilde->getMultiVector(0, bDomainThyraMode);
318 RCP<MultiVector> xtilde2 = bxtilde->getMultiVector(1, bDomainThyraMode);
321 RCP<MultiVector> xhat = MultiVectorFactory::Build(rcpX->getMap(), rcpX->getNumVectors());
322 RCP<BlockedMultiVector> bxhat = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(xhat);
323 RCP<MultiVector> xhat1 = bxhat->getMultiVector(0, bDomainThyraMode);
324 RCP<MultiVector> xhat2 = bxhat->getMultiVector(1, bDomainThyraMode);
330 residual->update(one, *rcpB, zero);
331 if(InitialGuessIsZero ==
false || run > 0)
332 A_->apply(*rcpX, *residual, Teuchos::NO_TRANS, -one, one);
336 xtilde1->putScalar(zero);
337 xtilde2->putScalar(zero);
338 velPredictSmoo_->Apply(*xtilde1, *r1);
342 RCP<MultiVector> schurCompRHS = rangeMapExtractor_->getVector(1, rcpB->getNumVectors(), bRangeThyraMode);
343 D_->apply(*xtilde1, *schurCompRHS);
345 schurCompRHS->update(one, *r2, -one);
349 schurCompSmoo_->Apply(*xtilde2, *schurCompRHS);
353 xhat2->update(omega, *xtilde2, zero);
356 RCP<MultiVector> xhat1_temp = domainMapExtractor_->getVector(0, rcpX->getNumVectors(), bDomainThyraMode);
357 G_->apply(*xhat2, *xhat1_temp);
359 xhat1->elementWiseMultiply(one, *diagFinv_, *xhat1_temp, zero);
360 xhat1->update(one, *xtilde1, -one);
362 rcpX->update(one, *bxhat, one);
365 if (bCopyResultX ==
true) {
366 RCP<const MultiVector> Xmerged = bX->Merge();
367 X.update(one, *Xmerged, zero);
372 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
373 RCP<MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
378 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
380 std::ostringstream out;
382 out <<
"{type = " << type_ <<
"}";
386 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
391 out0 <<
"Prec. type: " << type_ << std::endl;
394 if (verbLevel &
Debug) {
399 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
402 return Teuchos::OrdinalTraits<size_t>::invalid();
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
#define MUELU_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultScalar Scalar
virtual std::string description() const
Return a simple one-line description of this object.
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
virtual const Teuchos::ParameterList & GetParameterList() const =0
An exception safe way to call the method 'Level::SetFactoryManager()'.
SIMPLE smoother for 2x2 block matrices.
void SetVelocityPredictionFactoryManager(RCP< FactoryManager > FactManager)
void Setup(Level ¤tLevel)
Setup routine.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
std::string description() const
Return a simple one-line description of this object.
void DeclareInput(Level ¤tLevel) const
Input.
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
RCP< const ParameterList > GetValidParameterList() const
Input.
virtual ~SimpleSmoother()
Destructor.
RCP< SmootherPrototype > Copy() const
void AddFactoryManager(RCP< const FactoryManagerBase > FactManager, int pos)
Add a factory manager at a specific position.
void Apply(MultiVector &X, MultiVector const &B, bool InitialGuessIsZero=false) const
Apply the Braess Sarazin smoother.
void SetSchurCompFactoryManager(RCP< FactoryManager > FactManager)
SimpleSmoother()
Constructor.
bool IsSetup() const
Get the state of a smoother prototype.
static Teuchos::RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetLumpedMatrixDiagonal(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > const &A, const bool doReciprocal=false, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar tolReplacement=Teuchos::ScalarTraits< Scalar >::zero(), const bool replaceSingleEntryRowWithZero=false, const bool useAverageAbsDiagVal=false)
static Teuchos::RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetInverse(Teuchos::RCP< const Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > v, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar tolReplacement=Teuchos::ScalarTraits< Scalar >::zero())
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Debug
Print additional debugging information.
@ Parameters0
Print class parameters.