MueLu_RefMaxwell_decl.hpp

Go to the documentation of this file.

// @HEADER
//
// ***********************************************************************
//
//        MueLu: A package for multigrid based preconditioning
//                  Copyright 2012 Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact
//                    Jonathan Hu       (jhu@sandia.gov)
//                    Andrey Prokopenko (aprokop@sandia.gov)
//                    Ray Tuminaro      (rstumin@sandia.gov)
//
// ***********************************************************************
//
// @HEADER
#ifndef MUELU_REFMAXWELL_DECL_HPP
#define MUELU_REFMAXWELL_DECL_HPP

#include "MueLu_ConfigDefs.hpp"
#include "MueLu_BaseClass.hpp"

#include "MueLu_CoalesceDropFactory_fwd.hpp"
#include "MueLu_CoarseMapFactory_fwd.hpp"
#include "MueLu_CoordinatesTransferFactory_fwd.hpp"
#include "MueLu_TentativePFactory_fwd.hpp"
#include "MueLu_SaPFactory_fwd.hpp"
#include "MueLu_UncoupledAggregationFactory_fwd.hpp"
#include "MueLu_AggregationExportFactory_fwd.hpp"
#include "MueLu_Utilities_fwd.hpp"

#include "MueLu_MasterList.hpp"
#include "MueLu_FactoryManager.hpp"
#include "MueLu_SmootherBase_fwd.hpp"
#include "MueLu_SmootherPrototype_fwd.hpp"
#include "MueLu_PerfUtils_fwd.hpp"
#include "MueLu_RAPFactory_fwd.hpp"

#include "MueLu_CoalesceDropFactory_kokkos_fwd.hpp"
#include "MueLu_TentativePFactory_kokkos_fwd.hpp"

#include "MueLu_ZoltanInterface_fwd.hpp"
#include "MueLu_Zoltan2Interface_fwd.hpp"
#include "MueLu_RepartitionHeuristicFactory_fwd.hpp"
#include "MueLu_RepartitionFactory_fwd.hpp"
#include "MueLu_RebalanceAcFactory_fwd.hpp"
#include "MueLu_RebalanceTransferFactory_fwd.hpp"

#include "MueLu_SmootherFactory_fwd.hpp"
#include "MueLu_TrilinosSmoother_fwd.hpp"
#include "MueLu_Hierarchy_fwd.hpp"

#include "Xpetra_Operator.hpp"
#include "Xpetra_Map_fwd.hpp"
#include "Xpetra_Matrix_fwd.hpp"
#include "Xpetra_MapFactory_fwd.hpp"
#include "Xpetra_MatrixFactory_fwd.hpp"
#include "Xpetra_MultiVectorFactory_fwd.hpp"
#include "Xpetra_VectorFactory_fwd.hpp"
#include "Xpetra_CrsMatrixWrap_fwd.hpp"

namespace MueLu {

template <class Scalar,
          class LocalOrdinal,
          class GlobalOrdinal,
          class Node>
class RefMaxwell : public VerboseObject, public Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> {
#undef MUELU_REFMAXWELL_SHORT
#include "MueLu_UseShortNames.hpp"

 public:
  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitudeType;
  typedef typename Teuchos::ScalarTraits<Scalar>::coordinateType coordinateType;
  typedef typename Xpetra::MultiVector<coordinateType, LO, GO, NO> RealValuedMultiVector;

  RefMaxwell()
    : HierarchyCoarse11_(Teuchos::null)
    , Hierarchy22_(Teuchos::null)
    , disable_addon_(MasterList::getDefault<bool>("refmaxwell: disable addon"))
    , mode_(MasterList::getDefault<std::string>("refmaxwell: mode")) {
  }

  RefMaxwell(Teuchos::RCP<Hierarchy> HH, Teuchos::RCP<Hierarchy> H22)
    : HierarchyCoarse11_(HH)
    , Hierarchy22_(H22)
    , disable_addon_(MasterList::getDefault<bool>("refmaxwell: disable addon"))
    , mode_(MasterList::getDefault<std::string>("refmaxwell: mode")) {
  }

  RefMaxwell(const Teuchos::RCP<Matrix> &SM_Matrix,
             const Teuchos::RCP<Matrix> &Dk_1,
             const Teuchos::RCP<Matrix> &Dk_2,
             const Teuchos::RCP<Matrix> &D0,
             const Teuchos::RCP<Matrix> &M1_beta,
             const Teuchos::RCP<Matrix> &M1_alpha,
             const Teuchos::RCP<Matrix> &Mk_one,
             const Teuchos::RCP<Matrix> &Mk_1_one,
             const Teuchos::RCP<Matrix> &invMk_1_invBeta,
             const Teuchos::RCP<Matrix> &invMk_2_invAlpha,
             const Teuchos::RCP<MultiVector> &Nullspace11,
             const Teuchos::RCP<MultiVector> &Nullspace22,
             const Teuchos::RCP<RealValuedMultiVector> &NodalCoords,
             Teuchos::ParameterList &List,
             bool ComputePrec = true) {
    int spaceNumber = List.get<int>("refmaxwell: space number", 1);
    initialize(spaceNumber,
               Dk_1, Dk_2, D0,
               M1_beta, M1_alpha,
               Mk_one, Mk_1_one,
               invMk_1_invBeta, invMk_2_invAlpha,
               Nullspace11, Nullspace22, NodalCoords,
               List);
    resetMatrix(SM_Matrix, ComputePrec);
  }

  RefMaxwell(const Teuchos::RCP<Matrix> &SM_Matrix,
             const Teuchos::RCP<Matrix> &D0_Matrix,
             const Teuchos::RCP<Matrix> &Ms_Matrix,
             const Teuchos::RCP<Matrix> &M0inv_Matrix,
             const Teuchos::RCP<Matrix> &M1_Matrix,
             const Teuchos::RCP<MultiVector> &Nullspace11,
             const Teuchos::RCP<RealValuedMultiVector> &NodalCoords,
             Teuchos::ParameterList &List,
             bool ComputePrec = true) {
    initialize(D0_Matrix, Ms_Matrix, M0inv_Matrix, M1_Matrix, Nullspace11, NodalCoords, List);
    resetMatrix(SM_Matrix, ComputePrec);
  }

  RefMaxwell(const Teuchos::RCP<Matrix> &SM_Matrix,
             const Teuchos::RCP<Matrix> &D0_Matrix,
             const Teuchos::RCP<Matrix> &M0inv_Matrix,
             const Teuchos::RCP<Matrix> &M1_Matrix,
             const Teuchos::RCP<MultiVector> &Nullspace11,
             const Teuchos::RCP<RealValuedMultiVector> &NodalCoords,
             Teuchos::ParameterList &List,
             bool ComputePrec = true) {
    initialize(D0_Matrix, M1_Matrix, M0inv_Matrix, M1_Matrix, Nullspace11, NodalCoords, List);
    resetMatrix(SM_Matrix, ComputePrec);
  }

  RefMaxwell(const Teuchos::RCP<Matrix> &D0_Matrix,
             const Teuchos::RCP<Matrix> &M0inv_Matrix,
             const Teuchos::RCP<Matrix> &M1_Matrix,
             const Teuchos::RCP<MultiVector> &Nullspace11,
             const Teuchos::RCP<RealValuedMultiVector> &NodalCoords,
             Teuchos::ParameterList &List)
    : SM_Matrix_(Teuchos::null) {
    initialize(D0_Matrix, M1_Matrix, M0inv_Matrix, M1_Matrix, Nullspace11, NodalCoords, List);
  }

  RefMaxwell(const Teuchos::RCP<Matrix> &SM_Matrix,
             const Teuchos::RCP<Matrix> &D0_Matrix,
             const Teuchos::RCP<Matrix> &M1_Matrix,
             const Teuchos::RCP<MultiVector> &Nullspace11,
             const Teuchos::RCP<RealValuedMultiVector> &NodalCoords,
             Teuchos::ParameterList &List,
             bool ComputePrec) {
    initialize(D0_Matrix, M1_Matrix, Teuchos::null, M1_Matrix, Nullspace11, NodalCoords, List);
    resetMatrix(SM_Matrix, ComputePrec);
  }

  RefMaxwell(const Teuchos::RCP<Matrix> &D0_Matrix,
             const Teuchos::RCP<Matrix> &M1_Matrix,
             const Teuchos::RCP<MultiVector> &Nullspace11,
             const Teuchos::RCP<RealValuedMultiVector> &NodalCoords,
             Teuchos::ParameterList &List)
    : SM_Matrix_(Teuchos::null) {
    initialize(D0_Matrix, M1_Matrix, Teuchos::null, M1_Matrix, Nullspace11, NodalCoords, List);
  }

  RefMaxwell(const Teuchos::RCP<Matrix> &SM_Matrix,
             Teuchos::ParameterList &List,
             bool ComputePrec = true);

  virtual ~RefMaxwell() {}

  const Teuchos::RCP<const Map> getDomainMap() const;

  const Teuchos::RCP<const Map> getRangeMap() const;

  const Teuchos::RCP<Matrix> &getJacobian() const {
    return SM_Matrix_;
  }

  void setParameters(Teuchos::ParameterList &list);

  void compute(bool reuse = false);

  void resetMatrix(Teuchos::RCP<Matrix> SM_Matrix_new, bool ComputePrec = true);

  void apply(const MultiVector &X, MultiVector &Y,
             Teuchos::ETransp mode = Teuchos::NO_TRANS,
             Scalar alpha          = Teuchos::ScalarTraits<Scalar>::one(),
             Scalar beta           = Teuchos::ScalarTraits<Scalar>::zero()) const;

  bool hasTransposeApply() const;

  void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel = Teuchos::VERB_HIGH) const;

  void residual(const MultiVector &X,
                const MultiVector &B,
                MultiVector &R) const {
    using STS = Teuchos::ScalarTraits<Scalar>;
    R.update(STS::one(), B, STS::zero());
    this->apply(X, R, Teuchos::NO_TRANS, -STS::one(), STS::one());
  }

 private:
  Teuchos::RCP<Teuchos::ParameterList> getValidParamterList();

  void initialize(const Teuchos::RCP<Matrix> &D0_Matrix,
                  const Teuchos::RCP<Matrix> &Ms_Matrix,
                  const Teuchos::RCP<Matrix> &M0inv_Matrix,
                  const Teuchos::RCP<Matrix> &M1_Matrix,
                  const Teuchos::RCP<MultiVector> &Nullspace11,
                  const Teuchos::RCP<RealValuedMultiVector> &NodalCoords,
                  Teuchos::ParameterList &List);

  void initialize(const int k,
                  const Teuchos::RCP<Matrix> &Dk_1,
                  const Teuchos::RCP<Matrix> &Dk_2,
                  const Teuchos::RCP<Matrix> &D0,
                  const Teuchos::RCP<Matrix> &M1_beta,
                  const Teuchos::RCP<Matrix> &M1_alpha,
                  const Teuchos::RCP<Matrix> &Mk_one,
                  const Teuchos::RCP<Matrix> &Mk_1_one,
                  const Teuchos::RCP<Matrix> &invMk_1_invBeta,
                  const Teuchos::RCP<Matrix> &invMk_2_invAlpha,
                  const Teuchos::RCP<MultiVector> &Nullspace11,
                  const Teuchos::RCP<MultiVector> &Nullspace22,
                  const Teuchos::RCP<RealValuedMultiVector> &NodalCoords,
                  Teuchos::ParameterList &List);

  void determineSubHierarchyCommSizes(bool &doRebalancing, int &rebalanceStriding, int &numProcsCoarseA11, int &numProcsA22);

  void buildCoarse11Matrix();

  void rebalanceCoarse11Matrix(const int rebalanceStriding, const int numProcsCoarseA11);

  void build22Matrix(const bool reuse, const bool doRebalancing, const int rebalanceStriding, const int numProcsA22);

 public:  // due to Cuda build errors otherwise
  RCP<MultiVector> buildNullspace(const int spaceNumber, const Kokkos::View<bool *, typename Node::device_type> &bcs, const bool applyBCs);

  Teuchos::RCP<Matrix> buildProjection(const int spaceNumber, const RCP<MultiVector> &EdgeNullspace) const;

  void buildNodalProlongator(const Teuchos::RCP<Matrix> &A_nodal,
                             Teuchos::RCP<Matrix> &P_nodal,
                             Teuchos::RCP<MultiVector> &Nullspace_nodal,
                             Teuchos::RCP<RealValuedMultiVector> &Coords_nodal) const;

  RCP<Matrix> buildVectorNodalProlongator(const Teuchos::RCP<Matrix> &P_nodal) const;

  void buildProlongator(const int spaceNumber,
                        const Teuchos::RCP<Matrix> &A_nodal_Matrix,
                        const RCP<MultiVector> &EdgeNullspace,
                        Teuchos::RCP<Matrix> &edgeProlongator,
                        Teuchos::RCP<MultiVector> &coarseEdgeNullspace,
                        Teuchos::RCP<RealValuedMultiVector> &coarseNodalCoords) const;

  RCP<Matrix> buildAddon(const int spaceNumber);

 private:
  void setupSubSolve(Teuchos::RCP<Hierarchy> &hierarchy,
                     Teuchos::RCP<Operator> &thyraPrecOp,
                     const Teuchos::RCP<Matrix> &A,
                     const Teuchos::RCP<MultiVector> &Nullspace,
                     const Teuchos::RCP<RealValuedMultiVector> &Coords,
                     Teuchos::ParameterList &params,
                     std::string &label,
                     const bool reuse,
                     const bool isSingular = false);

  void setFineLevelSmoother11();

  void applyInverseAdditive(const MultiVector &RHS, MultiVector &X) const;

  void solveH(const MultiVector &RHS, MultiVector &X) const;

  void solve22(const MultiVector &RHS, MultiVector &X) const;

  void allocateMemory(int numVectors) const;

  void dump(const RCP<Matrix> &A, std::string name) const;

  void dump(const RCP<MultiVector> &X, std::string name) const;

  void dumpCoords(const RCP<RealValuedMultiVector> &X, std::string name) const;

  void dump(const Teuchos::ArrayRCP<bool> &v, std::string name) const;

  void dump(const Kokkos::View<bool *, typename Node::device_type> &v, std::string name) const;

  Teuchos::RCP<Teuchos::TimeMonitor> getTimer(std::string name, RCP<const Teuchos::Comm<int> > comm = Teuchos::null) const;

  Teuchos::RCP<Hierarchy> HierarchyCoarse11_, Hierarchy22_;
  Teuchos::RCP<SmootherBase> PreSmoother11_, PostSmoother11_;
  Teuchos::RCP<SmootherPrototype> PreSmootherData11_, PostSmootherData11_;
  RCP<Operator> thyraPrecOpH_, thyraPrecOp22_;
  int spaceNumber_;
  std::string solverName_;
  size_t dim_;
  Teuchos::RCP<Matrix> SM_Matrix_;
  Teuchos::RCP<Matrix> Dk_1_, Dk_1_T_;
  Teuchos::RCP<Matrix> Dk_2_;
  Teuchos::RCP<Matrix> D0_;
  Teuchos::RCP<Matrix> invMk_1_invBeta_, invMk_2_invAlpha_;
  Teuchos::RCP<Matrix> Mk_one_, Mk_1_one_;
  Teuchos::RCP<Matrix> M1_beta_, M1_alpha_;
  Teuchos::RCP<Matrix> P11_, R11_;
  Teuchos::RCP<Matrix> P22_, R22_;
  Teuchos::RCP<Matrix> coarseA11_, A22_, coarseA22_;
  Teuchos::RCP<Matrix> Addon11_;
  Teuchos::RCP<Matrix> Addon22_;
  Teuchos::RCP<const Map> DorigDomainMap_;
  Teuchos::RCP<const Import> DorigImporter_;
  Kokkos::View<bool *, typename Node::device_type::memory_space> BCrows11_, BCcols22_, BCdomain22_;
  int globalNumberBoundaryUnknowns11_, globalNumberBoundaryUnknowns22_;
  Teuchos::RCP<MultiVector> Nullspace11_, Nullspace22_;
  Teuchos::RCP<RealValuedMultiVector> NodalCoords_, CoordsCoarse11_, Coords22_;
  Teuchos::RCP<MultiVector> NullspaceCoarse11_;
  Teuchos::RCP<MultiVector> CoarseNullspace22_;
  Teuchos::RCP<const Import> ImporterCoarse11_, Importer22_;
  bool Dk_1_T_R11_colMapsMatch_;
  bool onlyBoundary11_, onlyBoundary22_;
  Teuchos::ParameterList parameterList_, precList11_, precList22_;
  Teuchos::RCP<Teuchos::ParameterList> coarseA11_AP_reuse_data_, coarseA11_RAP_reuse_data_;
  Teuchos::RCP<Teuchos::ParameterList> A22_AP_reuse_data_, A22_RAP_reuse_data_;
  bool disable_addon_, disable_addon_22_, dump_matrices_, useKokkos_, use_as_preconditioner_, implicitTranspose_, fuseProlongationAndUpdate_, syncTimers_, enable_reuse_, skipFirst11Level_, skipFirst22Level_;
  bool applyBCsToAnodal_, applyBCsToCoarse11_, applyBCsTo22_;
  int numItersCoarse11_, numIters22_;
  std::string mode_;

  mutable Teuchos::RCP<MultiVector> P11res_, P11x_, P11resSubComm_, P11xSubComm_, DresIntermediate_, Dres_, DxIntermediate_, Dx_, DresSubComm_, DxSubComm_, residual_, P11resTmp_, DresTmp_, DTR11Tmp_;
};

}  // namespace MueLu

#define MUELU_REFMAXWELL_SHORT
#endif  // MUELU_REFMAXWELL_DECL_HPP