MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_NullspaceFactory_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 MUELU_NULLSPACEFACTORY_DEF_HPP
11 #define MUELU_NULLSPACEFACTORY_DEF_HPP
12 
14 
15 #include <Xpetra_Matrix.hpp>
16 #include <Xpetra_MultiVectorFactory.hpp>
17 
18 #include "MueLu_Level.hpp"
19 #include "MueLu_Monitor.hpp"
20 #include "MueLu_MasterList.hpp"
21 #include "Xpetra_Access.hpp"
22 
23 namespace MueLu {
24 
25 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
27  RCP<ParameterList> validParamList = rcp(new ParameterList());
28 
29 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
30  SET_VALID_ENTRY("nullspace: calculate rotations");
31 #undef SET_VALID_ENTRY
32  validParamList->set<std::string>("Fine level nullspace", "Nullspace", "Variable name which is used to store null space multi vector on the finest level (default=\"Nullspace\"). For block matrices also \"Nullspace1\" to \"Nullspace9\" are accepted to describe the null space vectors for the (i,i) block (i=1..9).");
33 
34  validParamList->set<RCP<const FactoryBase>>("A", Teuchos::null, "Generating factory of the fine level matrix (only needed if default null space is generated)");
35  validParamList->set<RCP<const FactoryBase>>("Nullspace", Teuchos::null, "Generating factory of the fine level null space");
36  validParamList->set<RCP<const FactoryBase>>("Coordinates", Teuchos::null, "Generating factory of the coordinates");
37 
38  // TODO not very elegant.
39 
40  // 1/20/2016: we could add a sublist (e.g. "Nullspaces" which is excluded from parameter validation)
41  validParamList->set<RCP<const FactoryBase>>("Nullspace1", Teuchos::null, "Generating factory of the fine level null space associated with the (1,1) block in your n x n block matrix.");
42  validParamList->set<RCP<const FactoryBase>>("Nullspace2", Teuchos::null, "Generating factory of the fine level null space associated with the (2,2) block in your n x n block matrix.");
43  validParamList->set<RCP<const FactoryBase>>("Nullspace3", Teuchos::null, "Generating factory of the fine level null space associated with the (3,3) block in your n x n block matrix.");
44  validParamList->set<RCP<const FactoryBase>>("Nullspace4", Teuchos::null, "Generating factory of the fine level null space associated with the (4,4) block in your n x n block matrix.");
45  validParamList->set<RCP<const FactoryBase>>("Nullspace5", Teuchos::null, "Generating factory of the fine level null space associated with the (5,5) block in your n x n block matrix.");
46  validParamList->set<RCP<const FactoryBase>>("Nullspace6", Teuchos::null, "Generating factory of the fine level null space associated with the (6,6) block in your n x n block matrix.");
47  validParamList->set<RCP<const FactoryBase>>("Nullspace7", Teuchos::null, "Generating factory of the fine level null space associated with the (7,7) block in your n x n block matrix.");
48  validParamList->set<RCP<const FactoryBase>>("Nullspace8", Teuchos::null, "Generating factory of the fine level null space associated with the (8,8) block in your n x n block matrix.");
49  validParamList->set<RCP<const FactoryBase>>("Nullspace9", Teuchos::null, "Generating factory of the fine level null space associated with the (9,9) block in your n x n block matrix.");
50 
51  return validParamList;
52 }
53 
54 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
56  const ParameterList& pL = GetParameterList();
57  std::string nspName = pL.get<std::string>("Fine level nullspace");
58 
59  // only request "A" in DeclareInput if
60  // 1) there is not nspName (e.g. "Nullspace") is available in Level, AND
61  // 2) it is the finest level (i.e. LevelID == 0)
62  if (currentLevel.IsAvailable(nspName, NoFactory::get()) == false && currentLevel.GetLevelID() == 0)
63  Input(currentLevel, "A");
64 
65  if (currentLevel.GetLevelID() == 0 &&
66  currentLevel.IsAvailable("Coordinates", NoFactory::get()) && // we have coordinates (provided by user app)
67  pL.get<bool>("nullspace: calculate rotations")) { // and we want to calculate rotation modes
68  calculateRotations_ = true;
69  Input(currentLevel, "Coordinates");
70  }
71 
72  if (currentLevel.GetLevelID() != 0) {
73  // validate nullspaceFact_
74  // 1) The factory for "Nullspace" (or nspName) must not be Teuchos::null, since the default factory
75  // for "Nullspace" is a NullspaceFactory
76  // 2) The factory for "Nullspace" (or nspName) must be a TentativePFactory or any other TwoLevelFactoryBase derived object
77  // which generates the variable "Nullspace" as output
79  "MueLu::NullspaceFactory::DeclareInput(): You must declare an existing factory which "
80  "produces the variable \"Nullspace\" in the NullspaceFactory (e.g. a TentativePFactory).");
81  currentLevel.DeclareInput("Nullspace", GetFactory(nspName).get(), this); /* ! "Nullspace" and nspName mismatch possible here */
82  }
83 }
84 
85 template <class NullspaceType, class CoordsType, class MeanCoordsType, class LO>
87  private:
88  NullspaceType nullspace;
89  CoordsType coords;
90  MeanCoordsType mean;
93  typedef typename NullspaceType::value_type SC;
94 #if KOKKOS_VERSION >= 40799
95  typedef KokkosKernels::ArithTraits<SC> ATS;
96 #else
97  typedef Kokkos::ArithTraits<SC> ATS;
98 #endif
99 
100  public:
101  NullspaceFunctor(NullspaceType nullspace_, CoordsType coords_, MeanCoordsType mean_, LO numPDEs_, LO nullspaceDim_)
102  : nullspace(nullspace_)
103  , coords(coords_)
104  , mean(mean_)
105  , numPDEs(numPDEs_)
106  , nullspaceDim(nullspaceDim_) {
107  static_assert(static_cast<int>(NullspaceType::rank) == 2, "Nullspace needs to be a rank 2 view.");
108  static_assert(static_cast<int>(CoordsType::rank) == 2, "Coords needs to be a rank 2 view.");
109  static_assert(static_cast<int>(MeanCoordsType::rank) == 1, "Mean needs to be a rank 1 view.");
110  }
111 
112  KOKKOS_INLINE_FUNCTION
113  void operator()(const LO j) const {
114  SC one = ATS::one();
115  for (LO i = 0; i < numPDEs; i++)
116  nullspace(j * numPDEs + i, i) = one;
117  if ((nullspaceDim > numPDEs) && (numPDEs > 1)) {
118  // xy rotation
119  nullspace(j * numPDEs + 0, numPDEs) = -(coords(j, 1) - mean(1));
120  nullspace(j * numPDEs + 1, numPDEs) = (coords(j, 0) - mean(0));
121  }
122  if ((nullspaceDim == numPDEs + 3) && (numPDEs > 2)) {
123  // xz rotation
124  nullspace(j * numPDEs + 1, numPDEs + 1) = -(coords(j, 2) - mean(2));
125  nullspace(j * numPDEs + 2, numPDEs + 1) = (coords(j, 1) - mean(1));
126 
127  // yz rotation
128  nullspace(j * numPDEs + 0, numPDEs + 2) = -(coords(j, 2) - mean(2));
129  nullspace(j * numPDEs + 2, numPDEs + 2) = (coords(j, 0) - mean(0));
130  }
131  }
132 };
133 
134 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
136  FactoryMonitor m(*this, "Nullspace factory", currentLevel);
137 
138  RCP<MultiVector> nullspace;
139 
140  // TEUCHOS_TEST_FOR_EXCEPTION(currentLevel.GetLevelID() != 0, Exceptions::RuntimeError, "MueLu::NullspaceFactory::Build(): NullspaceFactory can be used for finest level (LevelID == 0) only.");
141  const ParameterList& pL = GetParameterList();
142  std::string nspName = pL.get<std::string>("Fine level nullspace");
143 
144  // get coordinates and compute mean of coordinates. (or centroid).
145 
147 
148  if (currentLevel.GetLevelID() == 0) {
149  if (currentLevel.IsAvailable(nspName, NoFactory::get())) {
150  // When a fine nullspace have already been defined by user using Set("Nullspace", ...) or
151  // Set("Nullspace1", ...), we use it.
152  nullspace = currentLevel.Get<RCP<MultiVector>>(nspName, NoFactory::get());
153  GetOStream(Runtime1) << "Use user-given nullspace " << nspName << ":"
154  << " nullspace dimension=" << nullspace->getNumVectors()
155  << " nullspace length=" << nullspace->getGlobalLength() << std::endl;
156 
157  } else {
158  // "Nullspace" (nspName) is not available
159  auto A = Get<RCP<Matrix>>(currentLevel, "A");
160 
161  // Determine numPDEs
162  LO numPDEs = 1;
163  if (A->IsView("stridedMaps") == true) {
164  Xpetra::viewLabel_t oldView = A->SwitchToView("stridedMaps"); // note: "stridedMaps are always non-overlapping (correspond to range and domain maps!)
165  TEUCHOS_TEST_FOR_EXCEPTION(rcp_dynamic_cast<const StridedMap>(A->getRowMap()).is_null(), Exceptions::BadCast,
166  "MueLu::CoalesceFactory::Build: cast to strided row map failed.");
167  numPDEs = rcp_dynamic_cast<const StridedMap>(A->getRowMap())->getFixedBlockSize();
168  oldView = A->SwitchToView(oldView);
169  }
170 
171  LO nullspaceDim = numPDEs;
172 
173  CoordsType coordsView;
174  MeanCoordsType meanView;
175  if (calculateRotations_) {
176  Coords = Get<RCP<RealValuedMultiVector>>(currentLevel, "Coordinates");
177  if (Coords->getNumVectors() > 1) nullspaceDim++;
178  if (Coords->getNumVectors() > 2) nullspaceDim += 2;
179 
180  meanView = MeanCoordsType("mean coords", Coords->getNumVectors());
181  Teuchos::Array<coordinate_type> hostMeans(Coords->getNumVectors());
182  Coords->meanValue(hostMeans);
183  Kokkos::View<typename RealValuedMultiVector::impl_scalar_type*, Kokkos::HostSpace> hostMeanView(hostMeans.getRawPtr(), hostMeans.size());
184  Kokkos::deep_copy(meanView, hostMeanView);
185  coordsView = Coords->getLocalViewDevice(Xpetra::Access::ReadOnly);
186  GetOStream(Runtime1) << "Generating nullspace with rotations: dimension = " << nullspaceDim << std::endl;
187  } else
188  GetOStream(Runtime1) << "Generating canonical nullspace: dimension = " << numPDEs << std::endl;
189 
190  nullspace = MultiVectorFactory::Build(A->getDomainMap(), nullspaceDim);
191 
192  fillNullspaceVector(nullspace, numPDEs, nullspaceDim, coordsView, meanView);
193  }
194 
195  } else {
196  // On coarser levels always use "Nullspace" as variable name, since it is expected by
197  // tentative P factory to be "Nullspace"
198  nullspace = currentLevel.Get<RCP<MultiVector>>("Nullspace", GetFactory(nspName).get()); // NOTE: "Nullspace" and nspName mismatch possible here
199  }
200 
201  // provide "Nullspace" variable on current level (used by TentativePFactory)
202  Set(currentLevel, "Nullspace", nullspace);
203 
204 } // Build
205 
206 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
208  RCP<BlockedMultiVector> bnsp = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(nullspace);
209  if (bnsp.is_null() == true) {
210  auto nullspaceView = nullspace->getLocalViewDevice(Xpetra::Access::OverwriteAll);
211 
212  int numBlocks = nullspace->getLocalLength() / numPDEs;
213  if (nullspaceDim > numPDEs)
214  TEUCHOS_TEST_FOR_EXCEPTION(numBlocks != coordsView.extent_int(0), Exceptions::RuntimeError, "MueLu::NullspaceFactory::fillNullspaceVector(): number of coordinates does not match ndofs/numPDEs.");
215 
216  NullspaceFunctor<decltype(nullspaceView), decltype(coordsView), decltype(meanView), LO> nullspaceFunctor(nullspaceView, coordsView, meanView, numPDEs, nullspaceDim);
217  Kokkos::parallel_for("MueLu:NullspaceF:Build:for", range_type(0, numBlocks), nullspaceFunctor);
218 
219  /*
220  // Scale columns to match what Galeri does. Not sure that this is necessary as the qr factorizatoin
221  // of the tentative prolongator also takes care of scaling issues. I'm leaving the code here
222  // just in case.
223  if ( (int) nullspaceDim > numPDEs ) {
224  Teuchos::Array<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> norms2(nullspaceDim);
225  nullspace->norm2(norms2);
226  Teuchos::Array<Scalar> norms2scalar(nullspaceDim);
227  for (int i = 0; i < nullspaceDim; i++)
228  norms2scalar[i] = norms2[0] / norms2[i];
229  nullspace->scale(norms2scalar);
230  }
231  */
232 
233  } else {
234  RCP<const BlockedMap> bmap = bnsp->getBlockedMap();
235  for (size_t r = 0; r < bmap->getNumMaps(); r++) {
236  Teuchos::RCP<MultiVector> part = bnsp->getMultiVector(r);
237  fillNullspaceVector(part, numPDEs, nullspaceDim, coordsView, meanView);
238  }
239  }
240 }
241 
242 } // namespace MueLu
243 
244 #endif // MUELU_NULLSPACEFACTORY_DEF_HPP
Exception indicating invalid cast attempted.
MueLu::DefaultLocalOrdinal LocalOrdinal
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level-&gt;Get&lt; RCP&lt;Matrix&gt; &gt;(&quot;A&quot;, factory) if factory == NULL =&gt; use default factory.
bool is_null(const boost::shared_ptr< T > &p)
RCP< const ParameterList > GetValidParameterList() const
Define valid parameters for internal factory parameters.
Kokkos::RangePolicy< local_ordinal_type, execution_space > range_type
T & get(const std::string &name, T def_value)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
typename RealValuedMultiVector::dual_view_type::t_dev_const_um CoordsType
LocalOrdinal LO
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
static const NoFactory * get()
void fillNullspaceVector(const RCP< MultiVector > &nullspace, LocalOrdinal numPDEs, LocalOrdinal nullspaceDim, CoordsType coordsView, MeanCoordsType meanView) const
void DeclareInput(Level &currentLevel) const
Specifies the data that this class needs, and the factories that generate that data.
Kokkos::View< typename RealValuedMultiVector::impl_scalar_type *, typename Node::memory_space > MeanCoordsType
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:63
KOKKOS_INLINE_FUNCTION void operator()(const LO j) const
std::string viewLabel_t
NullspaceType::value_type SC
NullspaceFunctor(NullspaceType nullspace_, CoordsType coords_, MeanCoordsType mean_, LO numPDEs_, LO nullspaceDim_)
Kokkos::ArithTraits< SC > ATS
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:51
Exception throws to report errors in the internal logical of the program.
Description of what is happening (more verbose)
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
#define SET_VALID_ENTRY(name)
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need&#39;s value has been saved.
void Build(Level &currentLevel) const
Build an object with this factory.
bool is_null() const