MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_UzawaSmoother_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#ifndef MUELU_UZAWASMOOTHER_DEF_HPP_
48#define MUELU_UZAWASMOOTHER_DEF_HPP_
49
50#include <Teuchos_ArrayViewDecl.hpp>
51#include <Teuchos_ScalarTraits.hpp>
52
53#include "MueLu_ConfigDefs.hpp"
54
55#include <Xpetra_Matrix.hpp>
56#include <Xpetra_BlockedCrsMatrix.hpp>
57#include <Xpetra_MultiVectorFactory.hpp>
58#include <Xpetra_VectorFactory.hpp>
59#include <Xpetra_ReorderedBlockedCrsMatrix.hpp>
60
62#include "MueLu_Level.hpp"
63#include "MueLu_Utilities.hpp"
64#include "MueLu_Monitor.hpp"
65#include "MueLu_HierarchyUtils.hpp"
67
68#include "MueLu_FactoryManager.hpp"
69
70namespace MueLu {
71
72 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
74 RCP<ParameterList> validParamList = rcp(new ParameterList());
75
76 validParamList->set<RCP<const FactoryBase>>("A", Teuchos::null, "Generating factory of the matrix A");
77 validParamList->set<Scalar> ("Damping factor", 1.0, "Damping/Scaling factor in SIMPLE");
78 validParamList->set<LocalOrdinal>("Sweeps", 1, "Number of SIMPLE sweeps (default = 1)");
79
80 return validParamList;
81 }
82
83 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
84 void UzawaSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::AddFactoryManager(RCP<const FactoryManagerBase> FactManager, int pos) {
85 TEUCHOS_TEST_FOR_EXCEPTION(pos < 0, Exceptions::RuntimeError, "MueLu::UzawaSmoother::AddFactoryManager: parameter \'pos\' must not be negative! error.");
86
87 size_t myPos = Teuchos::as<size_t>(pos);
88
89 if (myPos < FactManager_.size()) {
90 // replace existing entris in FactManager_ vector
91 FactManager_.at(myPos) = FactManager;
92 } else if( myPos == FactManager_.size()) {
93 // add new Factory manager in the end of the vector
94 FactManager_.push_back(FactManager);
95 } else { // if(myPos > FactManager_.size())
96 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
97 *out << "Warning: cannot add new FactoryManager at proper position " << pos << ". The FactoryManager is just appended to the end. Check this!" << std::endl;
98
99 // add new Factory manager in the end of the vector
100 FactManager_.push_back(FactManager);
101 }
102
103 }
104
105
106 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
108 AddFactoryManager(FactManager, 0); // overwrite factory manager for predicting the primary variable
109 }
110
111 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
113 AddFactoryManager(FactManager, 1); // overwrite factory manager for SchurComplement
114 }
115
116 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
118 currentLevel.DeclareInput("A",this->GetFactory("A").get());
119
120 TEUCHOS_TEST_FOR_EXCEPTION(FactManager_.size() != 2, Exceptions::RuntimeError,"MueLu::UzawaSmoother::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!");
121
122 // loop over all factory managers for the subblocks of blocked operator A
123 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
124 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
125 SetFactoryManager currentSFM (rcpFromRef(currentLevel), *it);
126
127 // request "Smoother" for current subblock row.
128 currentLevel.DeclareInput("PreSmoother",(*it)->GetFactory("Smoother").get());
129 }
130 }
131
132 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
134
135 FactoryMonitor m(*this, "Setup blocked Uzawa Smoother", currentLevel);
136
137 if (SmootherPrototype::IsSetup() == true)
138 this->GetOStream(Warnings0) << "MueLu::UzawaSmoother::Setup(): Setup() has already been called";
139
140 // extract blocked operator A from current level
141 A_ = Factory::Get<RCP<Matrix> > (currentLevel, "A");
142
143 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
144 TEUCHOS_TEST_FOR_EXCEPTION(bA == Teuchos::null, Exceptions::BadCast, "MueLu::UzawaSmoother::Setup: input matrix A is not of type BlockedCrsMatrix! error.");
145
146 // store map extractors
147 rangeMapExtractor_ = bA->getRangeMapExtractor();
148 domainMapExtractor_ = bA->getDomainMapExtractor();
149
150 // Store the blocks in local member variables
151 F_ = bA->getMatrix(0, 0);
152 G_ = bA->getMatrix(0, 1);
153 D_ = bA->getMatrix(1, 0);
154 Z_ = bA->getMatrix(1, 1);
155
156 // Set the Smoother
157 // carefully switch to the SubFactoryManagers (defined by the users)
158 {
159 RCP<const FactoryManagerBase> velpredictFactManager = FactManager_.at(0);
160 SetFactoryManager currentSFM (rcpFromRef(currentLevel), velpredictFactManager);
161 velPredictSmoo_ = currentLevel.Get< RCP<SmootherBase> >("PreSmoother",velpredictFactManager->GetFactory("Smoother").get());
162 }
163 {
164 RCP<const FactoryManagerBase> schurFactManager = FactManager_.at(1);
165 SetFactoryManager currentSFM (rcpFromRef(currentLevel), schurFactManager);
166 schurCompSmoo_ = currentLevel.Get< RCP<SmootherBase> >("PreSmoother", schurFactManager->GetFactory("Smoother").get());
167 }
168
170 }
171
172 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
173 void UzawaSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Apply(MultiVector& X, const MultiVector& B, bool /* InitialGuessIsZero */) const
174 {
175 TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError, "MueLu::UzawaSmoother::Apply(): Setup() has not been called");
176
177 Teuchos::RCP<Teuchos::FancyOStream> fos = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
178
179 SC zero = Teuchos::ScalarTraits<SC>::zero(), one = Teuchos::ScalarTraits<SC>::one();
180
181 bool bRangeThyraMode = rangeMapExtractor_->getThyraMode();
182 bool bDomainThyraMode = domainMapExtractor_->getThyraMode();
183
184 // extract parameters from internal parameter list
185 const ParameterList & pL = Factory::GetParameterList();
186 LocalOrdinal nSweeps = pL.get<LocalOrdinal>("Sweeps");
187 Scalar omega = pL.get<Scalar>("Damping factor");
188
189 // wrap current solution vector in RCP
190 RCP<MultiVector> rcpX = Teuchos::rcpFromRef(X);
191 RCP<const MultiVector> rcpB = Teuchos::rcpFromRef(B);
192
193 // make sure that both rcpX and rcpB are BlockedMultiVector objects
194 bool bCopyResultX = false;
195 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
196 MUELU_TEST_FOR_EXCEPTION(bA.is_null() == true, Exceptions::RuntimeError, "MueLu::BlockedGaussSeidelSmoother::Apply(): A_ must be a BlockedCrsMatrix");
197 RCP<BlockedMultiVector> bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
198 RCP<const BlockedMultiVector> bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
199
200 if(bX.is_null() == true) {
201 RCP<MultiVector> test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(),rcpX));
202 rcpX.swap(test);
203 bCopyResultX = true;
204 }
205
206 if(bB.is_null() == true) {
207 RCP<const MultiVector> test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB));
208 rcpB.swap(test);
209 }
210
211 // we now can guarantee that X and B are blocked multi vectors
212 bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
213 bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
214
215 // check the type of operator
216 RCP<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rbA = Teuchos::rcp_dynamic_cast<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(bA);
217 if(rbA.is_null() == false) {
218 // A is a ReorderedBlockedCrsMatrix
219 Teuchos::RCP<const Xpetra::BlockReorderManager > brm = rbA->getBlockReorderManager();
220
221 // check type of X vector
222 if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps()) {
223 // X is a blocked multi vector but incompatible to the reordered blocked operator A
224 Teuchos::RCP<MultiVector> test =
225 buildReorderedBlockedMultiVector(brm, bX);
226 rcpX.swap(test);
227 }
228 if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps()) {
229 // B is a blocked multi vector but incompatible to the reordered blocked operator A
230 Teuchos::RCP<const MultiVector> test =
231 buildReorderedBlockedMultiVector(brm, bB);
232 rcpB.swap(test);
233 }
234 }
235
236 // Throughout the rest of the algorithm rcpX and rcpB are used for solution vector and RHS
237
238 // create residual vector
239 // contains current residual of current solution X with rhs B
240 RCP<MultiVector> residual = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors());
241 RCP<BlockedMultiVector> bresidual = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(residual);
242 Teuchos::RCP<MultiVector> r1 = bresidual->getMultiVector(0,bRangeThyraMode);
243 Teuchos::RCP<MultiVector> r2 = bresidual->getMultiVector(1,bRangeThyraMode);
244
245 // helper vector 1
246 RCP<MultiVector> xtilde = MultiVectorFactory::Build(rcpX->getMap(), rcpX->getNumVectors());
247 RCP<BlockedMultiVector> bxtilde = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(xtilde);
248 RCP<MultiVector> xtilde1 = bxtilde->getMultiVector(0,bDomainThyraMode);
249 RCP<MultiVector> xtilde2 = bxtilde->getMultiVector(1,bDomainThyraMode);
250
251 // incrementally improve solution vector X
252 for (LocalOrdinal run = 0; run < nSweeps; ++run) {
253 // 1) calculate current residual
254 residual->update(one,*rcpB,zero); // residual = B
255 A_->apply(*rcpX, *residual, Teuchos::NO_TRANS, -one, one);
256
257 // 2) solve F * \Delta \tilde{x}_1 = r_1
258 // start with zero guess \Delta \tilde{x}_1
259 bxtilde->putScalar(zero);
260 velPredictSmoo_->Apply(*xtilde1,*r1);
261
262 // 3) calculate rhs for SchurComp equation
263 // r_2 - D \Delta \tilde{x}_1
264 RCP<MultiVector> schurCompRHS = MultiVectorFactory::Build(r2->getMap(),rcpB->getNumVectors());
265 D_->apply(*xtilde1,*schurCompRHS);
266 schurCompRHS->update(one,*r2,-one);
267
268 // 4) solve SchurComp equation
269 // start with zero guess \Delta \tilde{x}_2
270 schurCompSmoo_->Apply(*xtilde2,*schurCompRHS);
271
272 rcpX->update(omega,*bxtilde,one);
273 }
274
275 if (bCopyResultX == true) {
276 RCP<MultiVector> Xmerged = bX->Merge();
277 X.update(one, *Xmerged, zero);
278 }
279 }
280
281 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
282 RCP<MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
286
287 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
289 std::ostringstream out;
291 out << "{type = " << type_ << "}";
292 return out.str();
293 }
294
295 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
296 void UzawaSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::print(Teuchos::FancyOStream &out, const VerbLevel verbLevel) const {
298
299 if (verbLevel & Parameters0) {
300 out0 << "Prec. type: " << type_ << /*" Sweeps: " << nSweeps_ << " damping: " << omega_ <<*/ std::endl;
301 }
302
303 if (verbLevel & Debug) {
304 out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl;
305 }
306 }
307
308 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
310 // FIXME: This is a placeholder
311 return Teuchos::OrdinalTraits<size_t>::invalid();
312 }
313
314} // namespace MueLu
315
316
317#endif /* MUELU_UZAWASMOOTHER_DEF_HPP_ */
#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()'.
bool IsSetup() const
Get the state of a smoother prototype.
Block triangle Uzawa smoother for 2x2 block matrices.
RCP< SmootherPrototype > Copy() const
void Setup(Level &currentLevel)
Setup routine.
RCP< const ParameterList > GetValidParameterList() const
Input.
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.
std::string description() const
Return a simple one-line description of this object.
void SetSchurCompFactoryManager(RCP< FactoryManager > FactManager)
Set factory manager for internal SchurComplement handling.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
void DeclareInput(Level &currentLevel) const
Input.
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
void SetVelocityPredictionFactoryManager(RCP< FactoryManager > FactManager)
Set factory manager for internal velocity prediction.
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Debug
Print additional debugging information.
@ Parameters0
Print class parameters.