Teko_PCDStrategy.cpp

/*
// @HEADER
//
// ***********************************************************************
//
//      Teko: A package for block and physics based preconditioning
//                  Copyright 2010 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 Eric C. Cyr (eccyr@sandia.gov)
//
// ***********************************************************************
//
// @HEADER

*/

#include "Teko_PCDStrategy.hpp"

#include "Teuchos_TimeMonitor.hpp"
#include "Teko_Utilities.hpp"

namespace Teko {
namespace NS {

using Teuchos::TimeMonitor;

Teuchos::RCP<Teuchos::Time> PCDStrategy::initTimer_;
Teuchos::RCP<Teuchos::Time> PCDStrategy::invSTimer_;
Teuchos::RCP<Teuchos::Time> PCDStrategy::invFTimer_;
Teuchos::RCP<Teuchos::Time> PCDStrategy::opsTimer_;

void PCDStrategy::buildTimers() {
  if (initTimer_ == Teuchos::null)
    initTimer_ = TimeMonitor::getNewTimer("PCDStrategy::initializePrec");

  if (invSTimer_ == Teuchos::null)
    invSTimer_ = TimeMonitor::getNewTimer("PCDStrategy::initializePrec invS");

  if (invFTimer_ == Teuchos::null)
    invFTimer_ = TimeMonitor::getNewTimer("PCDStrategy::initializePrec invF");

  if (opsTimer_ == Teuchos::null)
    opsTimer_ = TimeMonitor::getNewTimer("PCDStrategy::initializePrec buildOps");
}

PCDStrategy::PCDStrategy() : massInverseType_(Diagonal), schurCompOrdering_(false) {
  pcdParams_ = Teuchos::rcp(new Teuchos::ParameterList);
  lapParams_ = Teuchos::rcp(new Teuchos::ParameterList);

  lapParams_->set("Name", getPressureLaplaceString());
  pcdParams_->set("Name", getPCDString());

  buildTimers();
}

PCDStrategy::PCDStrategy(const Teuchos::RCP<InverseFactory>& invFA,
                         const Teuchos::RCP<InverseFactory>& invS)
    : invFactoryF_(invFA),
      invFactoryS_(invS),
      massInverseType_(Diagonal),
      schurCompOrdering_(false) {
  pcdParams_ = Teuchos::rcp(new Teuchos::ParameterList);
  lapParams_ = Teuchos::rcp(new Teuchos::ParameterList);

  lapParams_->set("Name", getPressureLaplaceString());
  pcdParams_->set("Name", getPCDString());

  buildTimers();
}

const Teko::LinearOp PCDStrategy::getHatInvA00(const Teko::BlockedLinearOp& A,
                                               BlockPreconditionerState& state) const {
  initializeState(A, state);

  return state.getModifiableOp("invF");
}

const Teko::LinearOp PCDStrategy::getTildeInvA00(const Teko::BlockedLinearOp& A,
                                                 BlockPreconditionerState& state) const {
  initializeState(A, state);

  return state.getModifiableOp("invF");
}

const Teko::LinearOp PCDStrategy::getInvS(const Teko::BlockedLinearOp& A,
                                          BlockPreconditionerState& state) const {
  initializeState(A, state);

  return state.getLinearOp("invS");
}

void PCDStrategy::initializeState(const Teko::BlockedLinearOp& A,
                                  BlockPreconditionerState& state) const {
  Teko_DEBUG_SCOPE("PCDStrategy::initializeState", 10);
  TEUCHOS_ASSERT(getRequestHandler() != Teuchos::null);

  std::string pcdStr      = getPCDString();
  std::string presLapStr  = getPressureLaplaceString();
  std::string presMassStr = getPressureMassString();

  // no work to be done
  if (state.isInitialized()) return;

  Teuchos::TimeMonitor timer(*initTimer_, true);

  // extract sub blocks
  LinearOp F  = Teko::getBlock(0, 0, A);
  LinearOp Bt = Teko::getBlock(0, 1, A);
  LinearOp B  = Teko::getBlock(1, 0, A);
  LinearOp C  = Teko::getBlock(1, 1, A);

  LinearOp Qp = getRequestHandler()->request<LinearOp>(presMassStr);
  TEUCHOS_ASSERT(Qp != Teuchos::null);

  // build the inverse Laplacian complement
  LinearOp iQp;
  if (massInverseType_ == NotDiag) {
    ModifiableLinearOp& invMass = state.getModifiableOp("invMass");
    Teko_DEBUG_SCOPE("Building inv(Mass)", 10);

    if (invMass == Teuchos::null)
      invMass = buildInverse(*invFactoryS_, Qp);
    else
      rebuildInverse(*invFactoryS_, Qp, invMass);

    iQp = invMass;
  } else {
    Teko_DEBUG_MSG(
        "Building inverse mass of type \"" << Teko::getDiagonalName(massInverseType_) << "\"", 10);
    iQp = getInvDiagonalOp(Qp, massInverseType_);
  }

  // build the inverse Laplacian complement
  ModifiableLinearOp& invLaplace = state.getModifiableOp("invLaplace");
  {
    Teuchos::TimeMonitor timerInvS(*invSTimer_, true);

    // LinearOp laplace = getRequestHandler()->request<Teko::LinearOp>(presLapStr);
    LinearOp laplace = getRequestHandler()->request<Teko::LinearOp>(RequestMesg(lapParams_));
    TEUCHOS_ASSERT(laplace != Teuchos::null);
    if (invLaplace == Teuchos::null)
      invLaplace = buildInverse(*invFactoryS_, laplace);
    else
      rebuildInverse(*invFactoryS_, laplace, invLaplace);
  }

  // build the inverse Schur complement
  {
    Teko_DEBUG_SCOPE("Building S", 10);
    Teuchos::TimeMonitor timerS(*opsTimer_, true);

    // build Schur-complement
    // LinearOp pcd = getRequestHandler()->request<Teko::LinearOp>(pcdStr);
    LinearOp pcd = getRequestHandler()->request<Teko::LinearOp>(RequestMesg(pcdParams_));
    TEUCHOS_ASSERT(pcd != Teuchos::null);
    LinearOp invL = invLaplace;

    LinearOp invS;
    if (schurCompOrdering_ == false)
      invS = multiply(iQp, pcd, invL);
    else
      invS = multiply(invL, pcd, iQp);

    state.addLinearOp("invS", invS);
  }

  // build inverse F
  {
    Teko_DEBUG_SCOPE("Building inv(F)", 10);
    Teuchos::TimeMonitor timerInvF(*invFTimer_, true);

    ModifiableLinearOp& invF = state.getModifiableOp("invF");
    if (invF == Teuchos::null)
      invF = buildInverse(*invFactoryF_, F);
    else
      rebuildInverse(*invFactoryF_, F, invF);
  }

  // mark state as initialized
  state.setInitialized(true);
}

void PCDStrategy::initializeFromParameterList(const Teuchos::ParameterList& pl,
                                              const InverseLibrary& invLib) {
  Teko_DEBUG_SCOPE("PCDStrategy::initializeFromParameterList", 10);

  std::string invStr = "", invFStr = "", invSStr = "";
#if defined(Teko_ENABLE_Amesos)
  invStr = "Amesos";
#elif defined(Teko_ENABLE_Amesos2)
  invStr = "Amesos2";
#endif

  massInverseType_ = Diagonal;

  // "parse" the parameter list
  if (pl.isParameter("Inverse Type")) invStr = pl.get<std::string>("Inverse Type");
  if (pl.isParameter("Inverse F Type")) invFStr = pl.get<std::string>("Inverse F Type");
  if (pl.isParameter("Inverse Laplace Type")) invSStr = pl.get<std::string>("Inverse Laplace Type");
  if (pl.isParameter("Inverse Mass Type")) {
    std::string massInverseStr = pl.get<std::string>("Inverse Mass Type");

    // build inverse types
    massInverseType_ = getDiagonalType(massInverseStr);
  }
  if (pl.isParameter("Flip Schur Complement Ordering"))
    schurCompOrdering_ = pl.get<bool>("Flip Schur Complement Ordering");

  // set defaults as needed
  if (invFStr == "") invFStr = invStr;
  if (invSStr == "") invSStr = invStr;

  // read pressure laplace parameters
  if (pl.isSublist("Pressure Laplace Parameters"))
    lapParams_ =
        Teuchos::rcp(new Teuchos::ParameterList(pl.sublist("Pressure Laplace Parameters")));
  else
    lapParams_ = Teuchos::rcp(new Teuchos::ParameterList);

  // read pcd operator parameters
  if (pl.isSublist("Pressure Convection Diffusion Parameters"))
    pcdParams_ = Teuchos::rcp(
        new Teuchos::ParameterList(pl.sublist("Pressure Convection Diffusion Parameters")));
  else
    pcdParams_ = Teuchos::rcp(new Teuchos::ParameterList);

  // The user should not have already added this parameters
  TEUCHOS_TEST_FOR_EXCEPTION(
      lapParams_->isParameter("Name"), std::logic_error,
      "Teko: Parameter \"Name\" is not allowed in the sublist \"" + lapParams_->name() + "\"");
  TEUCHOS_TEST_FOR_EXCEPTION(
      lapParams_->isParameter("Tag"), std::logic_error,
      "Teko: Parameter \"Tag\" is not allowed in the sublist \"" + lapParams_->name() + "\"");
  TEUCHOS_TEST_FOR_EXCEPTION(
      pcdParams_->isParameter("Name"), std::logic_error,
      "Teko: Parameter \"Name\" is not allowed in the sublist \"" + pcdParams_->name() + "\"");
  TEUCHOS_TEST_FOR_EXCEPTION(
      pcdParams_->isParameter("Tag"), std::logic_error,
      "Teko: Parameter \"Tag\" is not allowed in the sublist \"" + pcdParams_->name() + "\"");

  Teko_DEBUG_MSG_BEGIN(5) DEBUG_STREAM << "PCD Strategy Parameters: " << std::endl;
  DEBUG_STREAM << "   inv type   = \"" << invStr << "\"" << std::endl;
  DEBUG_STREAM << "   inv F type = \"" << invFStr << "\"" << std::endl;
  DEBUG_STREAM << "   inv Laplace type = \"" << invSStr << "\"" << std::endl;
  DEBUG_STREAM << "   inv Mass type = \"" << Teko::getDiagonalName(massInverseType_) << "\""
               << std::endl;
  DEBUG_STREAM << "PCD Strategy Parameter list: " << std::endl;
  pl.print(DEBUG_STREAM);
  Teko_DEBUG_MSG_END()

      // build velocity inverse factory
      invFactoryF_ = invLib.getInverseFactory(invFStr);

  if (invFStr == invSStr)
    invFactoryS_ = invFactoryF_;
  else
    invFactoryS_ = invLib.getInverseFactory(invSStr);

  lapParams_->set("Name", getPressureLaplaceString());
  pcdParams_->set("Name", getPCDString());

  // setup a request for required operators
  getRequestHandler()->preRequest<Teko::LinearOp>(getPressureMassString());
  // getRequestHandler()->preRequest<Teko::LinearOp>(getPCDString());
  // getRequestHandler()->preRequest<Teko::LinearOp>(getPressureLaplaceString());
  getRequestHandler()->preRequest<Teko::LinearOp>(Teko::RequestMesg(lapParams_));
  getRequestHandler()->preRequest<Teko::LinearOp>(Teko::RequestMesg(pcdParams_));
}

Teuchos::RCP<Teuchos::ParameterList> PCDStrategy::getRequestedParameters() const {
  TEUCHOS_ASSERT(false);

  return Teuchos::null;
}

bool PCDStrategy::updateRequestedParameters(const Teuchos::ParameterList& /* pl */) {
  TEUCHOS_ASSERT(false);

  return true;
}

}  // end namespace NS
}  // end namespace Teko