MueLu Version of the Day
Loading...
Searching...
No Matches
Thyra_MueLuMaxwell1PreconditionerFactory_def.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// MueLu: A package for multigrid based preconditioning
4//
5// Copyright 2012 NTESS and the MueLu contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#ifndef THYRA_MUELU_MAXWELL1_PRECONDITIONER_FACTORY_DEF_HPP
11#define THYRA_MUELU_MAXWELL1_PRECONDITIONER_FACTORY_DEF_HPP
12
13#include <list>
14#include "Thyra_MueLuMaxwell1PreconditionerFactory_decl.hpp"
15#include <Xpetra_CrsMatrixWrap.hpp>
16#include <Xpetra_CrsMatrix.hpp>
17#include <Xpetra_Matrix.hpp>
18#include <Xpetra_ThyraUtils.hpp>
19#include <MueLu_Maxwell1.hpp>
20#include <Xpetra_TpetraHalfPrecisionOperator.hpp>
21
22#if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA)
23
24// This is not as general as possible, but should be good enough for most builds.
25#if ((defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) && !defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) && !defined(HAVE_TPETRA_INST_COMPLEX_FLOAT)) || \
26 (!defined(HAVE_TPETRA_INST_DOUBLE) && !defined(HAVE_TPETRA_INST_FLOAT) && defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) && defined(HAVE_TPETRA_INST_COMPLEX_FLOAT)) || \
27 (defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) && defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) && defined(HAVE_TPETRA_INST_COMPLEX_FLOAT)))
28#define MUELU_CAN_USE_MIXED_PRECISION
29#endif
30
31namespace Thyra {
32
33using Teuchos::ParameterList;
34using Teuchos::RCP;
35using Teuchos::rcp;
36using Teuchos::rcp_const_cast;
37using Teuchos::rcp_dynamic_cast;
38
39// Constructors/initializers/accessors
40
41template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
42MueLuMaxwell1PreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MueLuMaxwell1PreconditionerFactory()
43 : paramList_(rcp(new ParameterList())) {}
44
45// Overridden from PreconditionerFactoryBase
46
47template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
48bool MueLuMaxwell1PreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isCompatible(const LinearOpSourceBase<Scalar>& fwdOpSrc) const {
49 const RCP<const LinearOpBase<Scalar>> fwdOp = fwdOpSrc.getOp();
50
51 if (Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isTpetra(fwdOp)) return true;
52
53#ifdef HAVE_MUELU_EPETRA
54 if (Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isEpetra(fwdOp)) return true;
55#endif
56
57 return false;
58}
59
60template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
61RCP<PreconditionerBase<Scalar>> MueLuMaxwell1PreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::createPrec() const {
62 return Teuchos::rcp(new DefaultPreconditioner<Scalar>);
63}
64
65template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
66void MueLuMaxwell1PreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
67 initializePrec(const RCP<const LinearOpSourceBase<Scalar>>& fwdOpSrc, PreconditionerBase<Scalar>* prec, const ESupportSolveUse supportSolveUse) const {
68 // we are using typedefs here, since we are using objects from different packages (Xpetra, Thyra,...)
69 typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
70 typedef Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpThyUtils;
71 typedef Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpMat;
72 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
73 typedef Thyra::XpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyXpOp;
74#if defined(MUELU_CAN_USE_MIXED_PRECISION)
75 typedef Xpetra::TpetraHalfPrecisionOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpHalfPrecOp;
76 typedef Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpMV;
77 typedef typename XpHalfPrecOp::HalfScalar HalfScalar;
78 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
79 typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
80 typedef Xpetra::MultiVector<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> XphMV;
81 typedef Xpetra::MultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node> XpmMV;
82 typedef Xpetra::MultiVector<HalfMagnitude, LocalOrdinal, GlobalOrdinal, Node> XphmMV;
83 typedef Xpetra::Matrix<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> XphMat;
84#endif
85 Teuchos::TimeMonitor tM(*Teuchos::TimeMonitor::getNewTimer(std::string("ThyraMueLuMaxwell1::initializePrec")));
86
87 // Check precondition
88 TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
89 TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
90 TEUCHOS_ASSERT(prec);
91
92 // Create a copy, as we may remove some things from the list
93 ParameterList paramList = *paramList_;
94
95 // Retrieve wrapped concrete Xpetra matrix from FwdOp
96 const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
97 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
98
99 // Check whether it is Epetra/Tpetra
100 bool bIsEpetra = XpThyUtils::isEpetra(fwdOp);
101 bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
102 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == true && bIsTpetra == true));
103
104 // wrap the forward operator as an Xpetra::Matrix that MueLu can work with
105 // MueLu needs a non-const object as input
106 RCP<XpMat> A = XpThyUtils::toXpetra(Teuchos::rcp_const_cast<ThyLinOpBase>(fwdOp));
107 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A));
108
109 // Retrieve concrete preconditioner object
110 const Teuchos::Ptr<DefaultPreconditioner<Scalar>> defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar>*>(prec));
111 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
112
113 // extract preconditioner operator
114 RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
115 thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar>>(defaultPrec->getNonconstUnspecifiedPrecOp(), true);
116
117 // make a decision whether to (re)build the multigrid preconditioner or reuse the old one
118 // rebuild preconditioner if startingOver == true
119 // reuse preconditioner if startingOver == false
120 const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter("Maxwell1: enable reuse") || !paramList.get<bool>("Maxwell1: enable reuse"));
121 const bool useHalfPrecision = paramList.get<bool>("half precision", false) && bIsTpetra;
122
123 RCP<XpOp> xpPrecOp;
124 if (startingOver == true) {
125 // Convert to Xpetra
126 std::list<std::string> convertXpetra = {"Coordinates", "Nullspace", "Kn", "D0"};
127 for (auto it = convertXpetra.begin(); it != convertXpetra.end(); ++it)
128 Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(paramList, *it);
129
130 std::list<std::string> sublists = {"maxwell1: 11list", "maxwell1: 22list"};
131 for (auto itSublist = sublists.begin(); itSublist != sublists.end(); ++itSublist)
132 if (paramList.isSublist(*itSublist)) {
133 ParameterList& sublist = paramList.sublist(*itSublist);
134 for (int lvlNo = 0; lvlNo < 10; ++lvlNo) {
135 if (sublist.isSublist("level " + std::to_string(lvlNo) + " user data")) {
136 ParameterList& lvlList = sublist.sublist("level " + std::to_string(lvlNo) + " user data");
137 std::list<std::string> convertKeys;
138 for (auto it = lvlList.begin(); it != lvlList.end(); ++it)
139 convertKeys.push_back(lvlList.name(it));
140 for (auto it = convertKeys.begin(); it != convertKeys.end(); ++it)
141 Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(lvlList, *it);
142 }
143 }
144 }
145
146 ParameterList& sublist = paramList.sublist("maxwell1: 11list");
147 if (sublist.isParameter("D0")) {
148 Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(sublist, "D0");
149 }
150
151 paramList.set<bool>("Maxwell1: use as preconditioner", true);
152 if (useHalfPrecision) {
153#if defined(MUELU_CAN_USE_MIXED_PRECISION)
154
155 // convert to half precision
156 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
157 if (paramList.isType<RCP<XpmMV>>("Coordinates")) {
158 RCP<XpmMV> coords = paramList.get<RCP<XpmMV>>("Coordinates");
159 paramList.remove("Coordinates");
160 RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
161 paramList.set("Coordinates", halfCoords);
162 }
163 if (paramList.isType<RCP<XpMV>>("Nullspace")) {
164 RCP<XpMV> nullspace = paramList.get<RCP<XpMV>>("Nullspace");
165 paramList.remove("Nullspace");
166 RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
167 paramList.set("Nullspace", halfNullspace);
168 }
169 std::list<std::string> convertMat = {"Kn", "D0"};
170 for (auto it = convertMat.begin(); it != convertMat.end(); ++it) {
171 if (paramList.isType<RCP<XpMat>>(*it)) {
172 RCP<XpMat> M = paramList.get<RCP<XpMat>>(*it);
173 paramList.remove(*it);
174 RCP<XphMat> halfM = Xpetra::convertToHalfPrecision(M);
175 paramList.set(*it, halfM);
176 }
177 }
178
179 // build a new half-precision MueLu Maxwell1 preconditioner
180 RCP<MueLu::Maxwell1<HalfScalar, LocalOrdinal, GlobalOrdinal, Node>> halfPrec = rcp(new MueLu::Maxwell1<HalfScalar, LocalOrdinal, GlobalOrdinal, Node>(halfA, paramList, true));
181 xpPrecOp = rcp(new XpHalfPrecOp(halfPrec));
182#else
183 TEUCHOS_TEST_FOR_EXCEPT(true);
184#endif
185 } else {
186 // build a new MueLu Maxwell1 preconditioner
187 RCP<MueLu::Maxwell1<Scalar, LocalOrdinal, GlobalOrdinal, Node>> preconditioner = rcp(new MueLu::Maxwell1<Scalar, LocalOrdinal, GlobalOrdinal, Node>(A, paramList, true));
188 xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
189 }
190 } else {
191 // reuse old MueLu preconditioner stored in MueLu Xpetra operator and put in new matrix
192
193 RCP<ThyXpOp> thyXpOp = rcp_dynamic_cast<ThyXpOp>(thyra_precOp, true);
194 RCP<XpOp> xpOp = thyXpOp->getXpetraOperator();
195#if defined(MUELU_CAN_USE_MIXED_PRECISION)
196 RCP<XpHalfPrecOp> xpHalfPrecOp = rcp_dynamic_cast<XpHalfPrecOp>(xpOp);
197 if (!xpHalfPrecOp.is_null()) {
198 RCP<MueLu::Maxwell1<HalfScalar, LocalOrdinal, GlobalOrdinal, Node>> preconditioner = rcp_dynamic_cast<MueLu::Maxwell1<HalfScalar, LocalOrdinal, GlobalOrdinal, Node>>(xpHalfPrecOp->GetHalfPrecisionOperator(), true);
199 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
200 preconditioner->resetMatrix(halfA);
201 xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
202 } else
203#endif
204 {
205 RCP<MueLu::Maxwell1<Scalar, LocalOrdinal, GlobalOrdinal, Node>> preconditioner = rcp_dynamic_cast<MueLu::Maxwell1<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpOp, true);
206 preconditioner->resetMatrix(A);
207 xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
208 }
209 }
210
211 // wrap preconditioner in thyraPrecOp
212 RCP<const VectorSpaceBase<Scalar>> thyraRangeSpace = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpPrecOp->getRangeMap());
213 RCP<const VectorSpaceBase<Scalar>> thyraDomainSpace = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpPrecOp->getDomainMap());
214
215 RCP<ThyLinOpBase> thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace, xpPrecOp);
216 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraPrecOp));
217
218 defaultPrec->initializeUnspecified(thyraPrecOp);
219}
220
221template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
222void MueLuMaxwell1PreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
223 uninitializePrec(PreconditionerBase<Scalar>* prec, RCP<const LinearOpSourceBase<Scalar>>* fwdOp, ESupportSolveUse* supportSolveUse) const {
224 TEUCHOS_ASSERT(prec);
225
226 // Retrieve concrete preconditioner object
227 const Teuchos::Ptr<DefaultPreconditioner<Scalar>> defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar>*>(prec));
228 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
229
230 if (fwdOp) {
231 // TODO: Implement properly instead of returning default value
232 *fwdOp = Teuchos::null;
233 }
234
235 if (supportSolveUse) {
236 // TODO: Implement properly instead of returning default value
237 *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
238 }
239
240 defaultPrec->uninitialize();
241}
242
243// Overridden from ParameterListAcceptor
244template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
245void MueLuMaxwell1PreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setParameterList(RCP<ParameterList> const& paramList) {
246 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
247 paramList_ = paramList;
248}
249
250template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
251RCP<ParameterList> MueLuMaxwell1PreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getNonconstParameterList() {
252 return paramList_;
253}
254
255template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
256RCP<ParameterList> MueLuMaxwell1PreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::unsetParameterList() {
257 RCP<ParameterList> savedParamList = paramList_;
258 paramList_ = Teuchos::null;
259 return savedParamList;
260}
261
262template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
263RCP<const ParameterList> MueLuMaxwell1PreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getParameterList() const {
264 return paramList_;
265}
266
267template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
268RCP<const ParameterList> MueLuMaxwell1PreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getValidParameters() const {
269 static RCP<const ParameterList> validPL;
270
271 if (Teuchos::is_null(validPL))
272 validPL = rcp(new ParameterList());
273
274 return validPL;
275}
276
277// Public functions overridden from Teuchos::Describable
278template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
279std::string MueLuMaxwell1PreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::description() const {
280 return "Thyra::MueLuMaxwell1PreconditionerFactory";
281}
282} // namespace Thyra
283
284#endif // HAVE_MUELU_STRATIMIKOS
285
286#endif // ifdef THYRA_MUELU_MAXWELL1_PRECONDITIONER_FACTORY_DEF_HPP
Thyra_MueLuMaxwell1PreconditionerFactory_decl.hpp
MueLu::Maxwell1
Preconditioner (wrapped as a Xpetra::Operator) for Maxwell's equations in curl-curl form.
Definition MueLu_Maxwell1_decl.hpp:46
Thyra::XpetraLinearOp
Concrete Thyra::LinearOpBase subclass for Xpetra::Operator.
Definition Thyra_XpetraLinearOp_decl.hpp:28