MueLu_IndexManager_kokkos_def.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)
//                    Ray Tuminaro       (rstumin@sandia.gov)
//                    Luc Berger-Vergiat (lberge@sandia.gov)
//
// ***********************************************************************
//
// @HEADER
#ifndef MUELU_INDEXMANAGER_DEF_KOKKOS_HPP
#define MUELU_INDEXMANAGER_DEF_KOKKOS_HPP

#include <utility>

#include "Teuchos_OrdinalTraits.hpp"

#include <Xpetra_MapFactory.hpp>

#include "MueLu_ConfigDefs.hpp"
#include "MueLu_Types.hpp"
#include "MueLu_Exceptions.hpp"
#include <MueLu_IndexManager_kokkos_decl.hpp>

/*****************************************************************************

****************************************************************************/

namespace MueLu {

template <class LocalOrdinal, class GlobalOrdinal, class Node>
IndexManager_kokkos<LocalOrdinal, GlobalOrdinal, Node>::
    IndexManager_kokkos(const int NumDimensions,
                        const int interpolationOrder,
                        const int MyRank,
                        const ArrayView<const LO> LFineNodesPerDir,
                        const ArrayView<const int> CoarseRate)
  : myRank(MyRank)
  , coarseRate("coarsening rate")
  , endRate("endRate")
  , lFineNodesPerDir("lFineNodesPerDir")
  , coarseNodesPerDir("lFineNodesPerDir") {
  RCP<Teuchos::FancyOStream> out;
  if (const char* dbg = std::getenv("MUELU_INDEXMANAGER_DEBUG")) {
    out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
    out->setShowAllFrontMatter(false).setShowProcRank(true);
  } else {
    out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
  }

  setupIM(NumDimensions, interpolationOrder, CoarseRate, LFineNodesPerDir);

  *out << "Done setting up the IndexManager" << std::endl;

  computeMeshParameters();

  *out << "Computed Mesh Parameters" << std::endl;

}  // IndexManager_kokkos Constructor

template <class LocalOrdinal, class GlobalOrdinal, class Node>
void IndexManager_kokkos<LocalOrdinal, GlobalOrdinal, Node>::
    setupIM(const int NumDimensions, const int interpolationOrder,
            const ArrayView<const int> CoarseRate, const ArrayView<const LO> LFineNodesPerDir) {
  numDimensions       = NumDimensions;
  interpolationOrder_ = interpolationOrder;

  TEUCHOS_TEST_FOR_EXCEPTION((LFineNodesPerDir.size() != 3) && (LFineNodesPerDir.size() != numDimensions),
                             Exceptions::RuntimeError,
                             "LFineNodesPerDir has to be of size 3 or of size numDimensions!");

  typename Kokkos::View<LO[3], device_type>::HostMirror lFineNodesPerDir_h = Kokkos::create_mirror_view(lFineNodesPerDir);
  Kokkos::deep_copy(lFineNodesPerDir_h, lFineNodesPerDir);
  typename Kokkos::View<LO[3], device_type>::HostMirror coarseRate_h = Kokkos::create_mirror_view(coarseRate);
  Kokkos::deep_copy(coarseRate_h, coarseRate);

  // Load coarse rate, being careful about formating
  // Also load lFineNodesPerDir
  for (int dim = 0; dim < 3; ++dim) {
    if (dim < getNumDimensions()) {
      lFineNodesPerDir_h(dim) = LFineNodesPerDir[dim];
      if (CoarseRate.size() == 1) {
        coarseRate_h(dim) = CoarseRate[0];
      } else if (CoarseRate.size() == getNumDimensions()) {
        coarseRate_h(dim) = CoarseRate[dim];
      }
    } else {
      lFineNodesPerDir_h(dim) = 1;
      coarseRate_h(dim)       = 1;
    }
  }

  Kokkos::deep_copy(lFineNodesPerDir, lFineNodesPerDir_h);
  Kokkos::deep_copy(coarseRate, coarseRate_h);

}  // setupIM

template <class LocalOrdinal, class GlobalOrdinal, class Node>
void IndexManager_kokkos<LocalOrdinal, GlobalOrdinal, Node>::computeMeshParameters() {
  RCP<Teuchos::FancyOStream> out;
  if (const char* dbg = std::getenv("MUELU_INDEXMANAGER_DEBUG")) {
    out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
    out->setShowAllFrontMatter(false).setShowProcRank(true);
  } else {
    out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
  }

  typename Kokkos::View<int[3], device_type>::HostMirror coarseRate_h = Kokkos::create_mirror_view(coarseRate);
  typename Kokkos::View<int[3], device_type>::HostMirror endRate_h    = Kokkos::create_mirror_view(endRate);

  typename Kokkos::View<LO[3], device_type>::HostMirror lFineNodesPerDir_h  = Kokkos::create_mirror_view(lFineNodesPerDir);
  typename Kokkos::View<LO[3], device_type>::HostMirror coarseNodesPerDir_h = Kokkos::create_mirror_view(coarseNodesPerDir);
  Kokkos::deep_copy(lFineNodesPerDir_h, lFineNodesPerDir);
  Kokkos::deep_copy(coarseRate_h, coarseRate);

  lNumFineNodes10 = lFineNodesPerDir_h(1) * lFineNodesPerDir_h(0);
  lNumFineNodes   = lFineNodesPerDir_h(2) * lNumFineNodes10;
  for (int dim = 0; dim < 3; ++dim) {
    if (dim < numDimensions) {
      endRate_h(dim) = (lFineNodesPerDir_h(dim) - 1) % coarseRate_h(dim);
      if (endRate_h(dim) == 0) {
        endRate_h(dim) = coarseRate_h(dim);
      }

    } else {  // Default value for dim >= numDimensions
      endRate_h(dim) = 1;
    }
  }

  *out << "lFineNodesPerDir: {" << lFineNodesPerDir_h(0) << ", " << lFineNodesPerDir_h(1) << ", "
       << lFineNodesPerDir_h(2) << "}" << std::endl;
  *out << "endRate: {" << endRate_h(0) << ", " << endRate_h(1) << ", "
       << endRate_h(2) << "}" << std::endl;

  // Here one element can represent either the degenerate case of one node or the more general
  // case of two nodes, i.e. x---x is a 1D element with two nodes and x is a 1D element with
  // one node. This helps generating a 3D space from tensorial products...
  // A good way to handle this would be to generalize the algorithm to take into account the
  // discretization order used in each direction, at least in the FEM sense, since a 0 degree
  // discretization will have a unique node per element. This way 1D discretization can be
  // viewed as a 3D problem with one 0 degree element in the y direction and one 0 degre
  // element in the z direction.
  // !!! Operations below are aftecting both local and global values that have two         !!!
  // different orientations. Orientations can be interchanged using mapDirG2L and mapDirL2G.
  // coarseRate, endRate and offsets are in the global basis, as well as all the variables
  // starting with a g.
  // !!! while the variables starting with an l are in the local basis.                    !!!
  for (int dim = 0; dim < 3; ++dim) {
    if (dim < numDimensions) {
      // Check whether the partition includes the "end" of the mesh which means that endRate
      // will apply. Also make sure that endRate is not 0 which means that the mesh does not
      // require a particular treatment at the boundaries.
      coarseNodesPerDir_h(dim) = (lFineNodesPerDir_h(dim) - endRate_h(dim) - 1) / coarseRate_h(dim) + 2;

    } else {  // Default value for dim >= numDimensions
      // endRate[dim] = 1;
      coarseNodesPerDir_h(dim) = 1;
    }  // if (dim < numDimensions)

    // This would happen if the rank does not own any nodes but in that case a subcommunicator
    // should be used so this should really not be a concern.
    if (lFineNodesPerDir_h(dim) < 1) {
      coarseNodesPerDir_h(dim) = 0;
    }
  }  // Loop for dim=0:3

  // Compute cummulative values
  numCoarseNodes10 = coarseNodesPerDir_h(0) * coarseNodesPerDir_h(1);
  numCoarseNodes   = numCoarseNodes10 * coarseNodesPerDir_h(2);

  *out << "coarseNodesPerDir: {" << coarseNodesPerDir_h(0) << ", "
       << coarseNodesPerDir_h(1) << ", " << coarseNodesPerDir_h(2) << "}" << std::endl;
  *out << "numCoarseNodes=" << numCoarseNodes << std::endl;

  // Copy Host data to Device.
  Kokkos::deep_copy(coarseRate, coarseRate_h);
  Kokkos::deep_copy(endRate, endRate_h);
  Kokkos::deep_copy(lFineNodesPerDir, lFineNodesPerDir_h);
  Kokkos::deep_copy(coarseNodesPerDir, coarseNodesPerDir_h);
}

template <class LocalOrdinal, class GlobalOrdinal, class Node>
Array<LocalOrdinal> IndexManager_kokkos<LocalOrdinal, GlobalOrdinal, Node>::
    getCoarseNodesPerDirArray() const {
  typename LOTupleView::HostMirror coarseNodesPerDir_h = Kokkos::create_mirror_view(coarseNodesPerDir);
  Kokkos::deep_copy(coarseNodesPerDir_h, coarseNodesPerDir);
  Array<LO> coarseNodesPerDirArray(3);

  for (int dim = 0; dim < 3; ++dim) {
    coarseNodesPerDirArray[dim] = coarseNodesPerDir_h(dim);
  }

  return coarseNodesPerDirArray;
}  // getCoarseNodesData

}  // namespace MueLu

#define MUELU_INDEXMANAGER_KOKKOS_SHORT
#endif  // MUELU_INDEXMANAGER_DEF_KOKKOS_HPP