53 #ifndef MUELU_BRAESSSARAZINSMOOTHER_DEF_HPP_ 54 #define MUELU_BRAESSSARAZINSMOOTHER_DEF_HPP_ 56 #include "Teuchos_ArrayViewDecl.hpp" 57 #include "Teuchos_ScalarTraits.hpp" 61 #include <Xpetra_Matrix.hpp> 62 #include <Xpetra_CrsMatrixWrap.hpp> 63 #include <Xpetra_BlockedCrsMatrix.hpp> 64 #include <Xpetra_MultiVectorFactory.hpp> 65 #include <Xpetra_VectorFactory.hpp> 69 #include "MueLu_Utilities.hpp" 71 #include "MueLu_HierarchyHelpers.hpp" 75 #include "MueLu_SchurComplementFactory.hpp" 76 #include "MueLu_DirectSolver.hpp" 77 #include "MueLu_SmootherFactory.hpp" 78 #include "MueLu_FactoryManager.hpp" 82 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
84 : type_(
"Braess Sarazin"), A_(null)
94 SchurFact->SetParameter(
"omega",ParameterEntry(omega));
95 SchurFact->SetFactory(
"A", this->
GetFactory(
"A"));
98 ParameterList SCparams;
100 RCP<SmootherPrototype> smoProtoSC = rcp(
new DirectSolver(SCtype,SCparams) );
101 smoProtoSC->SetFactory(
"A", SchurFact);
103 RCP<SmootherFactory> SmooSCFact = rcp(
new SmootherFactory(smoProtoSC) );
106 FactManager->SetFactory(
"A", SchurFact);
107 FactManager->SetFactory(
"Smoother", SmooSCFact);
108 FactManager->SetIgnoreUserData(
true);
114 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
118 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
120 TEUCHOS_TEST_FOR_EXCEPTION(pos != 0,
Exceptions::RuntimeError,
"MueLu::BraessSarazinSmoother::AddFactoryManager: parameter \'pos\' must be zero! error.");
124 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
126 RCP<ParameterList> validParamList = rcp(
new ParameterList());
128 SC one = Teuchos::ScalarTraits<SC>::one();
130 validParamList->set<RCP<const FactoryBase> >(
"A", null,
"Generating factory of the matrix A");
132 validParamList->set<
bool> (
"lumping",
false,
"Use lumping to construct diag(A(0,0)), i.e. use row sum of the abs values on the diagonal " 133 "as approximation of A00 (and A00^{-1})");
134 validParamList->set<SC> (
"Damping factor", one,
"Damping/Scaling factor in BraessSarazin (usually has to be chosen > 1, default = 1.0)");
135 validParamList->set<LO> (
"Sweeps", 1,
"Number of BraessSarazin sweeps (default = 1)");
136 validParamList->set<ParameterList> (
"block1", ParameterList(),
"Sublist for parameters for SchurComplement block. At least has to contain some information about a smoother \"Smoother\" for variable \"A\" which is generated by a SchurComplementFactory.");
138 return validParamList;
141 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
143 this->
Input(currentLevel,
"A");
146 "MueLu::BraessSarazinSmoother::DeclareInput: FactManager_ must not be null! " 147 "Introduce a FactoryManager for the SchurComplement equation.");
166 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
171 this->
GetOStream(
Warnings0) <<
"MueLu::BreaessSarazinSmoother::Setup(): Setup() has already been called";
174 A_ = Factory::Get<RCP<Matrix> > (currentLevel,
"A");
175 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(
A_);
177 "MueLu::BraessSarazinSmoother::Setup: input matrix A is not of type BlockedCrsMatrix! error.");
190 A00_->CreateView(
"stridedMaps", bA->getRangeMap(0), bA->getDomainMap(0));
191 A01_->CreateView(
"stridedMaps", bA->getRangeMap(0), bA->getDomainMap(1));
192 A10_->CreateView(
"stridedMaps", bA->getRangeMap(1), bA->getDomainMap(0));
194 A11_->CreateView(
"stridedMaps", bA->getRangeMap(1), bA->getDomainMap(1));
197 SC omega = pL.get<SC>(
"Damping factor");
200 D_ = VectorFactory::Build(
A00_->getRowMap());
203 if (pL.get<
bool>(
"lumping") ==
false)
208 SC one = Teuchos::ScalarTraits<SC>::one();
210 ArrayRCP<SC> Ddata = D_->getDataNonConst(0);
211 for (GO row = 0; row < Ddata.size(); row++)
212 Ddata[row] = one / (diag[row]*omega);
218 smoo_ = currentLevel.Get<RCP<SmootherBase> >(
"PreSmoother",
FactManager_->GetFactory(
"Smoother").get());
219 S_ = currentLevel.Get<RCP<Matrix> > (
"A",
FactManager_->GetFactory(
"A").get());
225 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
228 "MueLu::BraessSarazinSmoother::Apply(): Setup() has not been called");
230 RCP<MultiVector> rcpX = rcpFromRef(X);
231 RCP<MultiVector> deltaX0 = MultiVectorFactory::Build(
A00_->getRowMap(), 1);
232 RCP<MultiVector> deltaX1 = MultiVectorFactory::Build(
A10_->getRowMap(), 1);
233 RCP<MultiVector> Rtmp = MultiVectorFactory::Build(
A10_->getRowMap(), 1);
235 typedef Teuchos::ScalarTraits<SC> STS;
236 SC one = STS::one(), zero = STS::zero();
240 LO nSweeps = pL.get<LO>(
"Sweeps");
243 if (InitialGuessIsZero) {
244 R = MultiVectorFactory::Build(B.getMap(), B.getNumVectors());
245 R->update(one, B, zero);
250 for (LO run = 0; run < nSweeps; ++run) {
259 deltaX0->putScalar(zero);
260 deltaX0->elementWiseMultiply(one, *
D_, *R0, zero);
261 A10_->apply(*deltaX0, *Rtmp);
262 Rtmp->update(one, *R1, -one);
266 deltaX1->putScalar(zero);
267 smoo_->Apply(*deltaX1, *Rtmp);
270 deltaX0->putScalar(zero);
271 A01_->apply(*deltaX1, *deltaX0);
272 deltaX0->update(one, *R0, -one);
274 deltaX0->elementWiseMultiply(one, *
D_, *R0, zero);
277 X0->update(one, *deltaX0, one);
278 X1->update(one, *deltaX1, one);
288 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
289 RCP<MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
294 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
297 std::ostringstream out;
299 out <<
"{type = " <<
type_ <<
"}";
303 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
308 out0 <<
"Prec. type: " <<
type_ << std::endl;
311 if (verbLevel &
Debug) {
static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
Important warning messages (one line)
Generic Smoother Factory for generating the smoothers of the MG hierarchy.
Exception indicating invalid cast attempted.
This class specifies the default factory that should generate some data on a Level if the data does n...
RCP< Matrix > A00_
matrices
BraessSarazin smoother for 2x2 block matrices.
RCP< Vector > D_
Inverse to approximation to block (0,0). Here, D_ = omega*inv(diag(A(0,0)))
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
RCP< const MapExtractorClass > domainMapExtractor_
domain map extractor (from A_ generated by AFact)
std::string type_
smoother type
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
RCP< Matrix > A_
block operator
Timer to be used in factories. Similar to Monitor but with additional timers.
Print additional debugging information.
static Teuchos::ArrayRCP< Scalar > GetLumpedMatrixDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
RCP< SmootherPrototype > Copy() const
Namespace for MueLu classes and methods.
virtual const Teuchos::ParameterList & GetParameterList() const
RCP< Matrix > S_
Schur complement.
RCP< const MapExtractorClass > rangeMapExtractor_
range map extractor (from A_ generated by AFact)
Class that encapsulates direct solvers. Autoselection of AmesosSmoother or Amesos2Smoother according ...
RCP< Matrix > A01_
Block (0,1) [typically, pressure gradient operator].
void DeclareInput(Level ¤tLevel) const
Input.
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &RHS)
Class that holds all level-specific information.
virtual ~BraessSarazinSmoother()
Destructor.
bool IsSetup() const
Get the state of a smoother prototype.
RCP< Matrix > A11_
Block (1,1) [typically, pressure stabilization term or null block].
RCP< const ParameterList > GetValidParameterList() const
Input.
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
Teuchos::FancyOStream & GetOStream(MsgType type, int thisProcRankOnly=0) const
Get an output stream for outputting the input message type.
BraessSarazinSmoother()
Constructor.
Factory for building the Schur Complement for a 2x2 block matrix.
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Crs2Op(RCP< Xpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
void AddFactoryManager(RCP< const FactoryManagerBase > FactManager, int pos=0)
Add a factory manager for BraessSarazin internal SchurComplement handling.
void Setup(Level ¤tLevel)
Setup routine.
RCP< Matrix > A10_
Block (1,0) [typically, divergence operator].
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the Braess Sarazin smoother.
Exception throws to report errors in the internal logical of the program.
An exception safe way to call the method 'Level::SetFactoryManager()'.
std::string description() const
Return a simple one-line description of this object.
Teuchos::RCP< SmootherBase > smoo_
Smoother for SchurComplement equation.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
void Input(Level &level, const std::string &varName) const
virtual std::string description() const
Return a simple one-line description of this object.
RCP< const FactoryManagerBase > FactManager_
Factory manager for creating the Schur Complement.
const RCP< const FactoryBase > GetFactory(const std::string &varName) const
Default implementation of FactoryAcceptor::GetFactory()