Teko_SmootherPreconditionerFactory.cpp

// @HEADER
// *****************************************************************************
//      Teko: A package for block and physics based preconditioning
//
// Copyright 2010 NTESS and the Teko contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

// Teko includes
#include "Teko_SmootherPreconditionerFactory.hpp"

#include "Teko_PreconditionerInverseFactory.hpp"
#include "Thyra_MultiVectorStdOps.hpp"

namespace Teko {

// A small utility class to help distinguish smoothers
class SmootherRequestMesg : public RequestMesg {
 public:
  SmootherRequestMesg(unsigned int block)
      : RequestMesg("__smoother_request_message__"), block_(block) {}

  unsigned int getBlock() const { return block_; }

 private:
  unsigned int block_;
};


SmootherLinearOp::SmootherLinearOp(const LinearOp& A, const LinearOp& invM,
                                   unsigned int applications, bool useDestAsInitialGuess)
    : A_(A),
      invM_(invM),
      applications_(applications),
      initialGuessType_(Unspecified),
      requestMesg_(Teuchos::null) {
  // set initial guess type
  initialGuessType_ = useDestAsInitialGuess ? DestAsInitialGuess : NoInitialGuess;
}

SmootherLinearOp::SmootherLinearOp(const LinearOp& A, const LinearOp& invM,
                                   unsigned int applications, unsigned int block)
    : A_(A),
      invM_(invM),
      applications_(applications),
      initialGuessType_(RequestInitialGuess),
      requestMesg_(Teuchos::null) {
  requestMesg_ = Teuchos::rcp(new SmootherRequestMesg(block));
}

void SmootherLinearOp::implicitApply(const MultiVector& b, MultiVector& x, const double alpha,
                                     const double beta) const {
  using Teuchos::RCP;

  MultiVector residual = deepcopy(b);  // residual = b
  MultiVector scrap    = deepcopy(b);  // scrap = b
  MultiVector error;                   // location for initial guess

  // construct initial guess: required to assign starting point for destination
  // vector appropriately
  switch (initialGuessType_) {
    case RequestInitialGuess:
      // get initial guess from request handler
      error = deepcopy(getRequestHandler()->request<MultiVector>(*requestMesg_));
      Thyra::assign<double>(x.ptr(), *error);  // x = error (initial guess)
      break;
    case DestAsInitialGuess:
      error = deepcopy(x);  // error = x
      break;
    case NoInitialGuess:
      Thyra::assign<double>(x.ptr(), 0.0);  // x = 0
      error = deepcopy(x);                  // error = x
      break;
    case Unspecified:
    default: TEUCHOS_ASSERT(false);
  }

  for (unsigned int current = 0; current < applications_; ++current) {
    // compute current residual
    Teko::applyOp(A_, error, scrap);
    Teko::update(-1.0, scrap, 1.0, residual);  // residual = residual-A*error

    // compute appoximate correction using residual
    Thyra::assign(error.ptr(), 0.0);  // set error to zero
    Teko::applyOp(invM_, residual, error);

    // update solution with error
    Teko::update(1.0, error, 1.0, x);  // x = x+error
  }
}

void SmootherLinearOp::setRequestHandler(const Teuchos::RCP<RequestHandler>& rh) {
  Teko_DEBUG_SCOPE("SmootherLinearOp::setRequestHandler", 10);
  requestHandler_ = rh;
}

Teuchos::RCP<RequestHandler> SmootherLinearOp::getRequestHandler() const {
  Teko_DEBUG_SCOPE("SmootherLinearOp::getRequestHandler", 10);
  return requestHandler_;
}

LinearOp buildSmootherLinearOp(const LinearOp& A, const LinearOp& invM, unsigned int applications,
                               bool useDestAsInitialGuess) {
  return Teuchos::rcp(new SmootherLinearOp(A, invM, applications, useDestAsInitialGuess));
}

LinearOp buildSmootherLinearOp(const LinearOp& A, const LinearOp& invM, unsigned int applications,
                               unsigned int initialGuessBlock) {
  return Teuchos::rcp(new SmootherLinearOp(A, invM, applications, initialGuessBlock));
}


SmootherPreconditionerFactory::SmootherPreconditionerFactory()
    : sweepCount_(0),
      initialGuessType_(Unspecified),
      initialGuessBlock_(0),
      precFactory_(Teuchos::null) {}

LinearOp SmootherPreconditionerFactory::buildPreconditionerOperator(
    LinearOp& lo, PreconditionerState& state) const {
  TEUCHOS_TEST_FOR_EXCEPTION(
      precFactory_ == Teuchos::null, std::runtime_error,
      "ERROR: Teko::SmootherPreconditionerFactory::buildPreconditionerOperator requires that a "
          << "preconditioner factory has been set. Currently it is null!");

  // preconditions
  TEUCHOS_ASSERT(sweepCount_ > 0);
  TEUCHOS_ASSERT(initialGuessType_ != Unspecified);
  TEUCHOS_ASSERT(precFactory_ != Teuchos::null);

  // build user specified preconditioner
  ModifiableLinearOp& invM = state.getModifiableOp("prec");
  if (invM == Teuchos::null)
    invM = Teko::buildInverse(*precFactory_, lo);
  else
    Teko::rebuildInverse(*precFactory_, lo, invM);

  // conditional on initial guess type, build the smoother
  switch (initialGuessType_) {
    case RequestInitialGuess:
      return buildSmootherLinearOp(lo, invM, sweepCount_, initialGuessBlock_);
    case DestAsInitialGuess:
      return buildSmootherLinearOp(lo, invM, sweepCount_, true);  // use an initial guess
    case NoInitialGuess:
      return buildSmootherLinearOp(lo, invM, sweepCount_, false);  // no initial guess
    case Unspecified:
    default: TEUCHOS_ASSERT(false);
  }

  // should never get here
  TEUCHOS_ASSERT(false);
  return Teuchos::null;
}

void SmootherPreconditionerFactory::initializeFromParameterList(
    const Teuchos::ParameterList& settings) {
  // declare all strings used by this initialization routine

  const std::string str_sweepCount         = "Sweep Count";
  const std::string str_initialGuessBlock  = "Initial Guess Block";
  const std::string str_destAsInitialGuess = "Destination As Initial Guess";
  const std::string str_precType           = "Preconditioner Type";

  // default parameters

  initialGuessType_  = Unspecified;
  initialGuessBlock_ = 0;
  sweepCount_        = 0;
  precFactory_       = Teuchos::null;

  // get sweep counts

  if (settings.isParameter(str_sweepCount)) sweepCount_ = settings.get<int>(str_sweepCount);

  // get initial guess

  if (settings.isParameter(str_initialGuessBlock)) {
    initialGuessBlock_ = settings.get<int>(str_initialGuessBlock);
    initialGuessType_  = RequestInitialGuess;
  }

  if (settings.isParameter(str_destAsInitialGuess)) {
    bool useDest = settings.get<bool>(str_destAsInitialGuess);
    if (useDest) {
      TEUCHOS_TEST_FOR_EXCEPTION(initialGuessType_ != Unspecified, std::runtime_error,
                                 "Cannot set both \"" << str_initialGuessBlock << "\" and \""
                                                      << str_destAsInitialGuess << "\"");

      initialGuessType_ = DestAsInitialGuess;
    }
  }

  // safe to assume if the other values are not turned on there is no initial guess
  if (initialGuessType_ == Unspecified) initialGuessType_ = NoInitialGuess;

  // get preconditioner factory

  TEUCHOS_TEST_FOR_EXCEPTION(
      not settings.isParameter(str_precType), std::runtime_error,
      "Parameter \"" << str_precType << "\" is required by a Teko::SmootherPreconditionerFactory");

  // grab library and preconditioner name
  Teuchos::RCP<const InverseLibrary> il = getInverseLibrary();
  std::string precName                  = settings.get<std::string>(str_precType);

  // build preconditioner factory
  precFactory_ = il->getInverseFactory(precName);
  TEUCHOS_TEST_FOR_EXCEPTION(
      precFactory_ == Teuchos::null, std::runtime_error,
      "ERROR: \"" << str_precType << "\" = " << precName << " could not be found");

  // post conditions required to be satisfied

  TEUCHOS_ASSERT(sweepCount_ > 0);
  TEUCHOS_ASSERT(initialGuessType_ != Unspecified);
  TEUCHOS_ASSERT(precFactory_ != Teuchos::null);
}

}  // end namespace Teko