46#ifndef MUELU_MAXWELL_UTILS_DEF_HPP
47#define MUELU_MAXWELL_UTILS_DEF_HPP
53#include "Xpetra_Map.hpp"
54#include "Xpetra_CrsMatrixUtils.hpp"
55#include "Xpetra_MatrixUtils.hpp"
60#include "MueLu_ThresholdAFilterFactory.hpp"
61#include "MueLu_Utilities.hpp"
62#include "MueLu_RAPFactory.hpp"
64#include "MueLu_Utilities_kokkos.hpp"
70 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
72 RCP<Matrix> & D0_Matrix_,
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_,
80 Teuchos::ArrayRCP<bool> & BCrows_,
81 Teuchos::ArrayRCP<bool> & BCcols_,
82 Teuchos::ArrayRCP<bool> & BCdomain_,
83 bool & allEdgesBoundary_,
84 bool & allNodesBoundary_) {
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());
101 auto BCrowsKokkos=BCrowsKokkos_;
102 Kokkos::parallel_reduce(BCrowsKokkos_.size(), KOKKOS_LAMBDA (
int i,
int & sum) {
107 auto BCdomainKokkos = BCdomainKokkos_;
108 Kokkos::parallel_reduce(BCdomainKokkos_.size(), KOKKOS_LAMBDA (
int i,
int & sum) {
109 if (BCdomainKokkos(i))
119 BCcols_.resize(D0_Matrix_->getColMap()->getLocalNumElements());
120 BCdomain_.resize(D0_Matrix_->getDomainMap()->getLocalNumElements());
123 for (
auto it = BCrows_.begin(); it != BCrows_.end(); ++it)
126 for (
auto it = BCdomain_.begin(); it != BCdomain_.end(); ++it)
131 MueLu_sumAll(SM_Matrix_->getRowMap()->getComm(), BCedgesLocal, BCedges_);
132 MueLu_sumAll(SM_Matrix_->getRowMap()->getComm(), BCnodesLocal, BCnodes_);
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();
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_) {
147 bool defaultFilter =
false;
152 if (parameterList_.get<
bool>(
"refmaxwell: filter D0",
true) && D0_Matrix_->getLocalMaxNumRowEntries()>2) {
156 fineLevel.
Set(
"A",D0_Matrix_);
157 fineLevel.
setlib(D0_Matrix_->getDomainMap()->lib());
160 fineLevel.
Request(
"A",ThreshFact.get());
161 ThreshFact->Build(fineLevel);
162 D0_Matrix_ = fineLevel.
Get< RCP<Matrix> >(
"A",ThreshFact.get());
165 defaultFilter =
true;
168 if (!M1_Matrix_.is_null() && parameterList_.get<
bool>(
"refmaxwell: filter M1", defaultFilter)) {
169 RCP<Vector> diag = VectorFactory::Build(M1_Matrix_->getRowMap());
171 M1_Matrix_->getLocalDiagCopy(*diag);
177 fineLevel.
Set(
"A",M1_Matrix_);
178 fineLevel.
setlib(M1_Matrix_->getDomainMap()->lib());
180 fineLevel.
Request(
"A",ThreshFact.get());
181 ThreshFact->Build(fineLevel);
182 M1_Matrix_ = fineLevel.
Get< RCP<Matrix> >(
"A",ThreshFact.get());
185 if (!Ms_Matrix_.is_null() && parameterList_.get<
bool>(
"refmaxwell: filter Ms", defaultFilter)) {
186 RCP<Vector> diag = VectorFactory::Build(Ms_Matrix_->getRowMap());
188 Ms_Matrix_->getLocalDiagCopy(*diag);
194 fineLevel.
Set(
"A",Ms_Matrix_);
195 fineLevel.
setlib(Ms_Matrix_->getDomainMap()->lib());
197 fineLevel.
Request(
"A",ThreshFact.get());
198 ThreshFact->Build(fineLevel);
199 Ms_Matrix_ = fineLevel.
Get< RCP<Matrix> >(
"A",ThreshFact.get());
202 if (!SM_Matrix_.is_null() && parameterList_.get<
bool>(
"refmaxwell: filter SM", defaultFilter)) {
203 RCP<Vector> diag = VectorFactory::Build(SM_Matrix_->getRowMap());
205 SM_Matrix_->getLocalDiagCopy(*diag);
211 fineLevel.
Set(
"A",SM_Matrix_);
212 fineLevel.
setlib(SM_Matrix_->getDomainMap()->lib());
214 fineLevel.
Request(
"A",ThreshFact.get());
215 ThreshFact->Build(fineLevel);
216 SM_Matrix_ = fineLevel.
Get< RCP<Matrix> >(
"A",ThreshFact.get());
223 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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);
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 ¶ms, std::string & label) {
239 Level fineLevel, coarseLevel;
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);
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);
260 coarseLevel.
Request(
"A", rapFact.get());
262 return coarseLevel.
Get< RCP<Matrix> >(
"A", rapFact.get());
272#define MUELU_MAXWELL_UTILS_SHORT
#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 ¶meterList, 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 ¶ms, 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.