MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_Maxwell_Utils_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#ifndef MUELU_MAXWELL_UTILS_DEF_HPP
47#define MUELU_MAXWELL_UTILS_DEF_HPP
48
49#include <sstream>
50
51#include "MueLu_ConfigDefs.hpp"
52
53#include "Xpetra_Map.hpp"
54#include "Xpetra_CrsMatrixUtils.hpp"
55#include "Xpetra_MatrixUtils.hpp"
56
58
59#include "MueLu_Level.hpp"
60#include "MueLu_ThresholdAFilterFactory.hpp"
61#include "MueLu_Utilities.hpp"
62#include "MueLu_RAPFactory.hpp"
63
64#include "MueLu_Utilities_kokkos.hpp"
65
66
67namespace MueLu {
68
69
70 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
72 RCP<Matrix> & D0_Matrix_,
73 magnitudeType rowSumTol,
74 bool useKokkos_,
75 Kokkos::View<bool*, typename Node::device_type> & BCrowsKokkos_,
76 Kokkos::View<bool*, typename Node::device_type> & BCcolsKokkos_,
77 Kokkos::View<bool*, typename Node::device_type> & BCdomainKokkos_,
78 int & BCedges_,
79 int & BCnodes_,
80 Teuchos::ArrayRCP<bool> & BCrows_,
81 Teuchos::ArrayRCP<bool> & BCcols_,
82 Teuchos::ArrayRCP<bool> & BCdomain_,
83 bool & allEdgesBoundary_,
84 bool & allNodesBoundary_) {
85 // clean rows associated with boundary conditions
86 // Find rows with only 1 or 2 nonzero entries, record them in BCrows_.
87 // BCrows_[i] is true, iff i is a boundary row
88 // BCcols_[i] is true, iff i is a boundary column
89 int BCedgesLocal = 0;
90 int BCnodesLocal = 0;
91 if (useKokkos_) {
92 BCrowsKokkos_ = Utilities_kokkos::DetectDirichletRows(*SM_Matrix_,Teuchos::ScalarTraits<magnitudeType>::eps(),/*count_twos_as_dirichlet=*/true);
93
94 if (rowSumTol > 0.)
95 Utilities_kokkos::ApplyRowSumCriterion(*SM_Matrix_, rowSumTol, BCrowsKokkos_);
96
97 BCcolsKokkos_ = Kokkos::View<bool*,typename Node::device_type>(Kokkos::ViewAllocateWithoutInitializing("dirichletCols"), D0_Matrix_->getColMap()->getLocalNumElements());
98 BCdomainKokkos_ = Kokkos::View<bool*,typename Node::device_type>(Kokkos::ViewAllocateWithoutInitializing("dirichletDomains"), D0_Matrix_->getDomainMap()->getLocalNumElements());
99 Utilities_kokkos::DetectDirichletColsAndDomains(*D0_Matrix_,BCrowsKokkos_,BCcolsKokkos_,BCdomainKokkos_);
100
101 auto BCrowsKokkos=BCrowsKokkos_;
102 Kokkos::parallel_reduce(BCrowsKokkos_.size(), KOKKOS_LAMBDA (int i, int & sum) {
103 if (BCrowsKokkos(i))
104 ++sum;
105 }, BCedgesLocal );
106
107 auto BCdomainKokkos = BCdomainKokkos_;
108 Kokkos::parallel_reduce(BCdomainKokkos_.size(), KOKKOS_LAMBDA (int i, int & sum) {
109 if (BCdomainKokkos(i))
110 ++sum;
111 }, BCnodesLocal);
112 } else
113 {
114 BCrows_ = Teuchos::arcp_const_cast<bool>(Utilities::DetectDirichletRows(*SM_Matrix_,Teuchos::ScalarTraits<magnitudeType>::eps(),/*count_twos_as_dirichlet=*/true));
115
116 if (rowSumTol > 0.)
117 Utilities::ApplyRowSumCriterion(*SM_Matrix_, rowSumTol, BCrows_);
118
119 BCcols_.resize(D0_Matrix_->getColMap()->getLocalNumElements());
120 BCdomain_.resize(D0_Matrix_->getDomainMap()->getLocalNumElements());
121 Utilities::DetectDirichletColsAndDomains(*D0_Matrix_,BCrows_,BCcols_,BCdomain_);
122
123 for (auto it = BCrows_.begin(); it != BCrows_.end(); ++it)
124 if (*it)
125 BCedgesLocal += 1;
126 for (auto it = BCdomain_.begin(); it != BCdomain_.end(); ++it)
127 if (*it)
128 BCnodesLocal += 1;
129 }
130
131 MueLu_sumAll(SM_Matrix_->getRowMap()->getComm(), BCedgesLocal, BCedges_);
132 MueLu_sumAll(SM_Matrix_->getRowMap()->getComm(), BCnodesLocal, BCnodes_);
133
134
135 allEdgesBoundary_ = Teuchos::as<Xpetra::global_size_t>(BCedges_) >= D0_Matrix_->getRangeMap()->getGlobalNumElements();
136 allNodesBoundary_ = Teuchos::as<Xpetra::global_size_t>(BCnodes_) >= D0_Matrix_->getDomainMap()->getGlobalNumElements();
137 }
138
139
140 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
142 RCP<Matrix> & D0_Matrix_,
143 RCP<Matrix> & SM_Matrix_,
144 RCP<Matrix> & M1_Matrix_,
145 RCP<Matrix> & Ms_Matrix_) {
146
147 bool defaultFilter = false;
148
149 // Remove zero entries from D0 if necessary.
150 // In the construction of the prolongator we use the graph of the
151 // matrix, so zero entries mess it up.
152 if (parameterList_.get<bool>("refmaxwell: filter D0", true) && D0_Matrix_->getLocalMaxNumRowEntries()>2) {
153 Level fineLevel;
154 fineLevel.SetFactoryManager(null);
155 fineLevel.SetLevelID(0);
156 fineLevel.Set("A",D0_Matrix_);
157 fineLevel.setlib(D0_Matrix_->getDomainMap()->lib());
158 // We expect D0 to have entries +-1, so any threshold value will do.
159 RCP<ThresholdAFilterFactory> ThreshFact = rcp(new ThresholdAFilterFactory("A",1.0e-8,/*keepDiagonal=*/false,/*expectedNNZperRow=*/2));
160 fineLevel.Request("A",ThreshFact.get());
161 ThreshFact->Build(fineLevel);
162 D0_Matrix_ = fineLevel.Get< RCP<Matrix> >("A",ThreshFact.get());
163
164 // If D0 has too many zeros, maybe SM and M1 do as well.
165 defaultFilter = true;
166 }
167
168 if (!M1_Matrix_.is_null() && parameterList_.get<bool>("refmaxwell: filter M1", defaultFilter)) {
169 RCP<Vector> diag = VectorFactory::Build(M1_Matrix_->getRowMap());
170 // find a reasonable absolute value threshold
171 M1_Matrix_->getLocalDiagCopy(*diag);
172 magnitudeType threshold = 1.0e-8 * diag->normInf();
173
174 Level fineLevel;
175 fineLevel.SetFactoryManager(null);
176 fineLevel.SetLevelID(0);
177 fineLevel.Set("A",M1_Matrix_);
178 fineLevel.setlib(M1_Matrix_->getDomainMap()->lib());
179 RCP<ThresholdAFilterFactory> ThreshFact = rcp(new ThresholdAFilterFactory("A",threshold,/*keepDiagonal=*/true));
180 fineLevel.Request("A",ThreshFact.get());
181 ThreshFact->Build(fineLevel);
182 M1_Matrix_ = fineLevel.Get< RCP<Matrix> >("A",ThreshFact.get());
183 }
184
185 if (!Ms_Matrix_.is_null() && parameterList_.get<bool>("refmaxwell: filter Ms", defaultFilter)) {
186 RCP<Vector> diag = VectorFactory::Build(Ms_Matrix_->getRowMap());
187 // find a reasonable absolute value threshold
188 Ms_Matrix_->getLocalDiagCopy(*diag);
189 magnitudeType threshold = 1.0e-8 * diag->normInf();
190
191 Level fineLevel;
192 fineLevel.SetFactoryManager(null);
193 fineLevel.SetLevelID(0);
194 fineLevel.Set("A",Ms_Matrix_);
195 fineLevel.setlib(Ms_Matrix_->getDomainMap()->lib());
196 RCP<ThresholdAFilterFactory> ThreshFact = rcp(new ThresholdAFilterFactory("A",threshold,/*keepDiagonal=*/true));
197 fineLevel.Request("A",ThreshFact.get());
198 ThreshFact->Build(fineLevel);
199 Ms_Matrix_ = fineLevel.Get< RCP<Matrix> >("A",ThreshFact.get());
200 }
201
202 if (!SM_Matrix_.is_null() && parameterList_.get<bool>("refmaxwell: filter SM", defaultFilter)) {
203 RCP<Vector> diag = VectorFactory::Build(SM_Matrix_->getRowMap());
204 // find a reasonable absolute value threshold
205 SM_Matrix_->getLocalDiagCopy(*diag);
206 magnitudeType threshold = 1.0e-8 * diag->normInf();
207
208 Level fineLevel;
209 fineLevel.SetFactoryManager(null);
210 fineLevel.SetLevelID(0);
211 fineLevel.Set("A",SM_Matrix_);
212 fineLevel.setlib(SM_Matrix_->getDomainMap()->lib());
213 RCP<ThresholdAFilterFactory> ThreshFact = rcp(new ThresholdAFilterFactory("A",threshold,/*keepDiagonal=*/true));
214 fineLevel.Request("A",ThreshFact.get());
215 ThreshFact->Build(fineLevel);
216 SM_Matrix_ = fineLevel.Get< RCP<Matrix> >("A",ThreshFact.get());
217 }
218
219 }
220
221
222
223 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
225 setMatvecParams(Matrix& A, RCP<ParameterList> matvecParams) {
226 RCP<const Import> xpImporter = A.getCrsGraph()->getImporter();
227 if (!xpImporter.is_null())
228 xpImporter->setDistributorParameters(matvecParams);
229 RCP<const Export> xpExporter = A.getCrsGraph()->getExporter();
230 if (!xpExporter.is_null())
231 xpExporter->setDistributorParameters(matvecParams);
232 }
233
234
235 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
236 RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
238 PtAPWrapper(RCP<Matrix>& A,RCP<Matrix>& P, ParameterList &params, std::string & label) {
239 Level fineLevel, coarseLevel;
240 fineLevel.SetFactoryManager(null);
241 coarseLevel.SetFactoryManager(null);
242 coarseLevel.SetPreviousLevel(rcpFromRef(fineLevel));
243 fineLevel.SetLevelID(0);
244 coarseLevel.SetLevelID(1);
245 fineLevel.Set("A",A);
246 coarseLevel.Set("P",P);
247 coarseLevel.setlib(A->getDomainMap()->lib());
248 fineLevel.setlib(A->getDomainMap()->lib());
249 coarseLevel.setObjectLabel(label);
250 fineLevel.setObjectLabel(label);
251
252 RCP<RAPFactory> rapFact = rcp(new RAPFactory());
253 ParameterList rapList = *(rapFact->GetValidParameterList());
254 rapList.set("transpose: use implicit", true);
255 rapList.set("rap: fix zero diagonals", params.get<bool>("rap: fix zero diagonals", true));
256 rapList.set("rap: fix zero diagonals threshold", params.get<double>("rap: fix zero diagonals threshold", Teuchos::ScalarTraits<double>::eps()));
257 rapList.set("rap: triple product", params.get<bool>("rap: triple product", false));
258 rapFact->SetParameterList(rapList);
259
260 coarseLevel.Request("A", rapFact.get());
261
262 return coarseLevel.Get< RCP<Matrix> >("A", rapFact.get());
263 }
264
265
266
267
268
269
270} // namespace
271
272#define MUELU_MAXWELL_UTILS_SHORT
273#endif //ifdef MUELU_MAXWELL_UTILS_DEF_HPP
#define MueLu_sumAll(rcpComm, in, out)
Class that holds all level-specific information.
void setlib(Xpetra::UnderlyingLib lib2)
void SetLevelID(int levelID)
Set level number.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
void Request(const FactoryBase &factory)
Increment the storage counter for all the inputs of a factory.
void SetPreviousLevel(const RCP< Level > &previousLevel)
void SetFactoryManager(const RCP< const FactoryManagerBase > &factoryManager)
Set default factories (used internally by Hierarchy::SetLevel()).
static void removeExplicitZeros(Teuchos::ParameterList &parameterList, RCP< Matrix > &D0_Matrix, RCP< Matrix > &SM_Matrix, RCP< Matrix > &M1_Matrix, RCP< Matrix > &Ms_Matrix)
Remove explicit zeros.
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitudeType
static void detectBoundaryConditionsSM(RCP< Matrix > &SM_Matrix, RCP< Matrix > &D0_Matrix, magnitudeType rowSumTol, bool useKokkos_, Kokkos::View< bool *, typename Node::device_type > &BCrowsKokkos, Kokkos::View< bool *, typename Node::device_type > &BCcolsKokkos, Kokkos::View< bool *, typename Node::device_type > &BCdomainKokkos, int &BCedges, int &BCnodes, Teuchos::ArrayRCP< bool > &BCrows, Teuchos::ArrayRCP< bool > &BCcols, Teuchos::ArrayRCP< bool > &BCdomain, bool &allEdgesBoundary, bool &allNodesBoundary)
Detect Dirichlet boundary conditions.
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > PtAPWrapper(RCP< Matrix > &A, RCP< Matrix > &P, Teuchos::ParameterList &params, std::string &label)
static void setMatvecParams(Matrix &A, RCP< ParameterList > matvecParams)
Sets matvec params on a matrix.
Factory for building coarse matrices.
Factory for building a thresholded operator.
static void DetectDirichletColsAndDomains(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Kokkos::View< bool *, typename Node::device_type > &dirichletRows, Kokkos::View< bool *, typename Node::device_type > dirichletCols, Kokkos::View< bool *, typename Node::device_type > dirichletDomain)
Detects Dirichlet columns & domains from a list of Dirichlet rows.
static void ApplyRowSumCriterion(const Matrix &A, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType rowSumTol, Kokkos::View< bool *, typename NO::device_type > &dirichletRows)
static Kokkos::View< bool *, typename NO::device_type > DetectDirichletRows(const Matrix &A, const Magnitude &tol=Teuchos::ScalarTraits< typename Teuchos::ScalarTraits< SC >::magnitudeType >::zero(), const bool count_twos_as_dirichlet=false)
Detect Dirichlet rows.
static void ApplyRowSumCriterion(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude rowSumTol, Teuchos::ArrayRCP< bool > &dirichletRows)
static void DetectDirichletColsAndDomains(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Teuchos::ArrayRCP< bool > &dirichletRows, Teuchos::ArrayRCP< bool > dirichletCols, Teuchos::ArrayRCP< bool > dirichletDomains)
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::magnitude(0.), const bool count_twos_as_dirichlet=false)
Namespace for MueLu classes and methods.