MueLu  Version of the Day
MueLu_BraessSarazinSmoother_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 /*
47  * MueLu_BraessSarazinSmoother_def.hpp
48  *
49  * Created on: Apr 16, 2012
50  * Author: wiesner
51  */
52 
53 #ifndef MUELU_BRAESSSARAZINSMOOTHER_DEF_HPP_
54 #define MUELU_BRAESSSARAZINSMOOTHER_DEF_HPP_
55 
56 #include "Teuchos_ArrayViewDecl.hpp"
57 #include "Teuchos_ScalarTraits.hpp"
58 
59 #include "MueLu_ConfigDefs.hpp"
60 
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>
66 
68 #include "MueLu_Level.hpp"
69 #include "MueLu_Utilities.hpp"
70 #include "MueLu_Monitor.hpp"
71 #include "MueLu_HierarchyHelpers.hpp"
72 #include "MueLu_SmootherBase.hpp"
73 
74 // include files for default FactoryManager
75 #include "MueLu_SchurComplementFactory.hpp"
76 #include "MueLu_DirectSolver.hpp"
77 #include "MueLu_SmootherFactory.hpp"
78 #include "MueLu_FactoryManager.hpp"
79 
80 namespace MueLu {
81 
82  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
84  : type_("Braess Sarazin"), A_(null)
85  {
86  //Factory::SetParameter("Sweeps", ParameterEntry(sweeps));
87  //Factory::SetParameter("Damping factor",ParameterEntry(omega));
88 
89 #if 0
90  // when declaring default factories without overwriting them leads to a multipleCallCheck exception
91  // TODO: debug into this
92  // workaround: always define your factory managers outside either using the C++ API or the XML files
93  RCP<SchurComplementFactory> SchurFact = rcp(new SchurComplementFactory());
94  SchurFact->SetParameter("omega",ParameterEntry(omega));
95  SchurFact->SetFactory("A", this->GetFactory("A"));
96 
97  // define smoother/solver for BraessSarazin
98  ParameterList SCparams;
99  std::string SCtype;
100  RCP<SmootherPrototype> smoProtoSC = rcp( new DirectSolver(SCtype,SCparams) );
101  smoProtoSC->SetFactory("A", SchurFact);
102 
103  RCP<SmootherFactory> SmooSCFact = rcp( new SmootherFactory(smoProtoSC) );
104 
105  RCP<FactoryManager> FactManager = rcp(new FactoryManager());
106  FactManager->SetFactory("A", SchurFact);
107  FactManager->SetFactory("Smoother", SmooSCFact);
108  FactManager->SetIgnoreUserData(true);
109 
110  AddFactoryManager(FactManager,0);
111 #endif
112  }
113 
114  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
116 
118  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
119  void BraessSarazinSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::AddFactoryManager(RCP<const FactoryManagerBase> FactManager, int pos) {
120  TEUCHOS_TEST_FOR_EXCEPTION(pos != 0, Exceptions::RuntimeError, "MueLu::BraessSarazinSmoother::AddFactoryManager: parameter \'pos\' must be zero! error.");
121  FactManager_ = FactManager;
122  }
123 
124  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
126  RCP<ParameterList> validParamList = rcp(new ParameterList());
127 
128  SC one = Teuchos::ScalarTraits<SC>::one();
129 
130  validParamList->set<RCP<const FactoryBase> >("A", null, "Generating factory of the matrix A");
131 
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.");
137 
138  return validParamList;
139  }
140 
141  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
143  this->Input(currentLevel, "A");
144 
145  TEUCHOS_TEST_FOR_EXCEPTION(FactManager_.is_null(), Exceptions::RuntimeError,
146  "MueLu::BraessSarazinSmoother::DeclareInput: FactManager_ must not be null! "
147  "Introduce a FactoryManager for the SchurComplement equation.");
148 
149  // carefully call DeclareInput after switching to internal FactoryManager
150  {
151  SetFactoryManager currentSFM(rcpFromRef(currentLevel), FactManager_);
152 
153  // request "Smoother" for current subblock row.
154  currentLevel.DeclareInput("PreSmoother", FactManager_->GetFactory("Smoother").get());
155 
156  // request Schur matrix just in case
157  currentLevel.DeclareInput("A", FactManager_->GetFactory("A").get());
158  }
159  }
160 
161  // Setup routine can be summarized in 4 steps:
162  // - set the map extractors
163  // - set the blocks
164  // - create and set the inverse of the diagonal of F
165  // - set the smoother for the Schur Complement
166  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
168  FactoryMonitor m(*this, "Setup Smoother", currentLevel);
169 
170  if (SmootherPrototype::IsSetup() == true)
171  this->GetOStream(Warnings0) << "MueLu::BreaessSarazinSmoother::Setup(): Setup() has already been called";
172 
173  // Extract blocked operator A from current level
174  A_ = Factory::Get<RCP<Matrix> > (currentLevel, "A");
175  RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
176  TEUCHOS_TEST_FOR_EXCEPTION(bA.is_null(), Exceptions::BadCast,
177  "MueLu::BraessSarazinSmoother::Setup: input matrix A is not of type BlockedCrsMatrix! error.");
178 
179  // Store map extractors
180  rangeMapExtractor_ = bA->getRangeMapExtractor();
181  domainMapExtractor_ = bA->getDomainMapExtractor();
182 
183  // Store the blocks in local member variables
188 
189  // TODO move this to BlockedCrsMatrix->getMatrix routine...
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));
193  if (!A11_.is_null())
194  A11_->CreateView("stridedMaps", bA->getRangeMap(1), bA->getDomainMap(1));
195 
196  const ParameterList& pL = Factory::GetParameterList();
197  SC omega = pL.get<SC>("Damping factor");
198 
199  // Create the inverse of the diagonal of F
200  D_ = VectorFactory::Build(A00_->getRowMap());
201 
202  ArrayRCP<SC> diag;
203  if (pL.get<bool>("lumping") == false)
205  else
207 
208  SC one = Teuchos::ScalarTraits<SC>::one();
209 
210  ArrayRCP<SC> Ddata = D_->getDataNonConst(0);
211  for (GO row = 0; row < Ddata.size(); row++)
212  Ddata[row] = one / (diag[row]*omega);
213 
214  // Set the Smoother
215  // carefully switch to the SubFactoryManagers (defined by the users)
216  {
217  SetFactoryManager currentSFM(rcpFromRef(currentLevel), FactManager_);
218  smoo_ = currentLevel.Get<RCP<SmootherBase> >("PreSmoother", FactManager_->GetFactory("Smoother").get());
219  S_ = currentLevel.Get<RCP<Matrix> > ("A", FactManager_->GetFactory("A").get());
220  }
221 
223  }
224 
225  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
226  void BraessSarazinSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Apply(MultiVector& X, const MultiVector& B, bool InitialGuessIsZero) const {
227  TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError,
228  "MueLu::BraessSarazinSmoother::Apply(): Setup() has not been called");
229 
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);
234 
235  typedef Teuchos::ScalarTraits<SC> STS;
236  SC one = STS::one(), zero = STS::zero();
237 
238  // extract parameters from internal parameter list
239  const ParameterList& pL = Factory::GetParameterList();
240  LO nSweeps = pL.get<LO>("Sweeps");
241 
242  RCP<MultiVector> R;
243  if (InitialGuessIsZero) {
244  R = MultiVectorFactory::Build(B.getMap(), B.getNumVectors());
245  R->update(one, B, zero);
246  } else {
247  R = Utilities::Residual(*A_, X, B);
248  }
249 
250  for (LO run = 0; run < nSweeps; ++run) {
251  // Extract corresponding subvectors from X and R
252  RCP<MultiVector> R0 = rangeMapExtractor_ ->ExtractVector(R, 0);
253  RCP<MultiVector> R1 = rangeMapExtractor_ ->ExtractVector(R, 1);
254 
255  RCP<MultiVector> X0 = domainMapExtractor_->ExtractVector(rcpX, 0);
256  RCP<MultiVector> X1 = domainMapExtractor_->ExtractVector(rcpX, 1);
257 
258  // Calculate Rtmp = R1 - D * deltaX0 (equation 8.14)
259  deltaX0->putScalar(zero);
260  deltaX0->elementWiseMultiply(one, *D_, *R0, zero); // deltaX0 = D * R0 (equation 8.13)
261  A10_->apply(*deltaX0, *Rtmp); // Rtmp = A10*D*deltaX0 (intermediate step)
262  Rtmp->update(one, *R1, -one); // Rtmp = R1 - A10*D*deltaX0
263 
264  // Compute deltaX1 (pressure correction)
265  // We use user provided preconditioner
266  deltaX1->putScalar(zero); // just for safety
267  smoo_->Apply(*deltaX1, *Rtmp);
268 
269  // Compute deltaX0
270  deltaX0->putScalar(zero); // just for safety
271  A01_->apply(*deltaX1, *deltaX0); // deltaX0 = A01*deltaX1
272  deltaX0->update(one, *R0, -one); // deltaX0 = R0 - A01*deltaX1
273  R0.swap(deltaX0);
274  deltaX0->elementWiseMultiply(one, *D_, *R0, zero); // deltaX0 = D*(R0 - A01*deltaX1)
275 
276  // Update solution
277  X0->update(one, *deltaX0, one);
278  X1->update(one, *deltaX1, one);
279 
280  domainMapExtractor_->InsertVector(X0, 0, rcpX);
281  domainMapExtractor_->InsertVector(X1, 1, rcpX);
282 
283  if (run < nSweeps-1)
284  R = Utilities::Residual(*A_, X, B);
285  }
286  }
287 
288  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
289  RCP<MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
291  return rcp (new BraessSarazinSmoother (*this));
292  }
293 
294  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
296  description () const {
297  std::ostringstream out;
299  out << "{type = " << type_ << "}";
300  return out.str();
301  }
302 
303  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
304  void BraessSarazinSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::print(Teuchos::FancyOStream &out, const VerbLevel verbLevel) const {
306 
307  if (verbLevel & Parameters0) {
308  out0 << "Prec. type: " << type_ << /*" Sweeps: " << nSweeps_ << " damping: " << omega_ <<*/ std::endl;
309  }
310 
311  if (verbLevel & Debug) {
312  out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl;
313  }
314  }
315 
316 } // namespace MueLu
317 
318 #endif /* MUELU_BRAESSSARAZINSMOOTHER_DEF_HPP_ */
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...
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)
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
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< 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 &currentLevel) 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.
Definition: MueLu_Level.hpp:99
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.
Print class parameters.
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 &currentLevel)
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 &#39;Level::SetFactoryManager()&#39;.
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()