Teko_LSCPreconditionerFactory.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_LSCPreconditionerFactory.hpp"

#include "Thyra_DefaultMultipliedLinearOp.hpp"
#include "Thyra_DefaultAddedLinearOp.hpp"
#include "Thyra_DefaultIdentityLinearOp.hpp"
#include "Thyra_DefaultZeroLinearOp.hpp"

#include "Teko_LU2x2InverseOp.hpp"
#include "Teko_Utilities.hpp"
#include "Teko_BlockUpperTriInverseOp.hpp"
#include "Teko_StaticLSCStrategy.hpp"
#include "Teko_InvLSCStrategy.hpp"
#include "Teko_PresLaplaceLSCStrategy.hpp"
#include "Teko_LSCSIMPLECStrategy.hpp"

#include "Teuchos_Time.hpp"

namespace Teko {
namespace NS {

using Teuchos::rcp;
using Teuchos::RCP;
using Teuchos::rcp_dynamic_cast;

using Thyra::add;
using Thyra::identity;
using Thyra::multiply;

// Stabilized constructor
LSCPreconditionerFactory::LSCPreconditionerFactory(const LinearOp& invF, const LinearOp& invBQBtmC,
                                                   const LinearOp& invD, const LinearOp& invMass)
    : invOpsStrategy_(rcp(new StaticLSCStrategy(invF, invBQBtmC, invD, invMass))),
      isSymmetric_(true) {}

// Stable constructor
LSCPreconditionerFactory::LSCPreconditionerFactory(const LinearOp& invF, const LinearOp& invBQBtmC,
                                                   const LinearOp& invMass)
    : invOpsStrategy_(rcp(new StaticLSCStrategy(invF, invBQBtmC, invMass))), isSymmetric_(true) {}

// fully generic constructor
LSCPreconditionerFactory::LSCPreconditionerFactory(const RCP<LSCStrategy>& strategy)
    : invOpsStrategy_(strategy), isSymmetric_(true) {}

LSCPreconditionerFactory::LSCPreconditionerFactory() : isSymmetric_(true) {}

// for PreconditionerFactoryBase

// initialize a newly created preconditioner object
LinearOp LSCPreconditionerFactory::buildPreconditionerOperator(
    BlockedLinearOp& blockOp, BlockPreconditionerState& state) const {
  Teko_DEBUG_SCOPE("LSCPreconditionerFactory::buildPreconditionerOperator", 10);
  Teko_DEBUG_EXPR(Teuchos::Time timer(""));
  Teko_DEBUG_EXPR(Teuchos::Time totalTimer(""));
  Teko_DEBUG_EXPR(totalTimer.start());

  // extract sub-matrices from source operator
  LinearOp F  = blockOp->getBlock(0, 0);
  LinearOp B  = blockOp->getBlock(1, 0);
  LinearOp Bt = blockOp->getBlock(0, 1);

  if (not isSymmetric_) Bt = scale(-1.0, adjoint(B));

  // build what is neccessary for the state object
  Teko_DEBUG_EXPR(timer.start(true));
  invOpsStrategy_->buildState(blockOp, state);
  Teko_DEBUG_EXPR(timer.stop());
  Teko_DEBUG_MSG("LSCPrecFact::buildPO BuildStateTime = " << timer.totalElapsedTime(), 2);

  // extract operators from strategy
  Teko_DEBUG_EXPR(timer.start(true));
  LinearOp invF      = invOpsStrategy_->getInvF(blockOp, state);
  LinearOp invBQBtmC = invOpsStrategy_->getInvBQBt(blockOp, state);
  LinearOp invBHBtmC = invOpsStrategy_->getInvBHBt(blockOp, state);
  LinearOp outerStab = invOpsStrategy_->getOuterStabilization(blockOp, state);
  LinearOp innerStab = invOpsStrategy_->getInnerStabilization(blockOp, state);

  // if necessary build an identity mass matrix
  LinearOp invMass  = invOpsStrategy_->getInvMass(blockOp, state);
  LinearOp HScaling = invOpsStrategy_->getHScaling(blockOp, state);
  if (invMass == Teuchos::null) invMass = identity<double>(F->range());
  if (HScaling == Teuchos::null) HScaling = identity<double>(F->range());
  Teko_DEBUG_EXPR(timer.stop());
  Teko_DEBUG_MSG("LSCPrecFact::buildPO GetInvTime = " << timer.totalElapsedTime(), 2);

  // need to build Schur complement,  inv(P) = inv(B*Bt)*(B*F*Bt)*inv(B*Bt)

  // first construct middle operator: M = B * inv(Mass) * F * inv(Mass) * Bt
  LinearOp M =
      //          (B * inv(Mass) ) * F * (inv(Mass) * Bt)
      multiply(multiply(B, invMass), F, multiply(HScaling, Bt));
  if (innerStab != Teuchos::null)  // deal with inner stabilization
    M = add(M, innerStab);

  // now construct a linear operator schur complement
  LinearOp invPschur;
  if (outerStab != Teuchos::null)
    invPschur = add(multiply(invBQBtmC, M, invBHBtmC), outerStab);
  else
    invPschur = multiply(invBQBtmC, M, invBHBtmC);

  // build the preconditioner operator: Use LDU or upper triangular approximation
  if (invOpsStrategy_->useFullLDU()) {
    Teko_DEBUG_EXPR(totalTimer.stop());
    Teko_DEBUG_MSG("LSCPrecFact::buildPO TotalTime = " << totalTimer.totalElapsedTime(), 2);

    // solve using a full LDU decomposition
    return createLU2x2InverseOp(blockOp, invF, invPschur, "LSC-LDU");
  } else {
    // build diagonal operations
    std::vector<LinearOp> invDiag(2);
    invDiag[0] = invF;
    invDiag[1] = Thyra::scale(-1.0, invPschur);

    // get upper triangular matrix
    BlockedLinearOp U = getUpperTriBlocks(blockOp);

    Teko_DEBUG_EXPR(totalTimer.stop());
    Teko_DEBUG_MSG("LSCPrecFact::buildPO TotalTime = " << totalTimer.totalElapsedTime(), 2);

    // solve using only one inversion of F
    return createBlockUpperTriInverseOp(U, invDiag, "LSC-Upper");
  }
}

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

  RCP<const InverseLibrary> invLib = getInverseLibrary();

  if (pl.isParameter("Is Symmetric")) isSymmetric_ = pl.get<bool>("Is Symmetric");

  std::string name = "Basic Inverse";
  if (pl.isParameter("Strategy Name")) name = pl.get<std::string>("Strategy Name");
  const Teuchos::ParameterEntry* pe = pl.getEntryPtr("Strategy Settings");

  // check for a mistake in input file
  if (name != "Basic Inverse" && pe == 0) {
    RCP<Teuchos::FancyOStream> out = getOutputStream();
    *out << "LSC Construction failed: ";
    *out << "Strategy \"" << name << "\" requires a \"Strategy Settings\" sublist" << std::endl;
    throw std::runtime_error("LSC Construction failed: Strategy Settings not set");
  }

  // get the parameter list to construct the strategy
  Teuchos::RCP<const Teuchos::ParameterList> stratPL = Teuchos::rcpFromRef(pl);
  if (pe != 0) stratPL = Teuchos::rcpFromRef(pl.sublist("Strategy Settings"));

  // build the strategy object
  RCP<LSCStrategy> strategy = buildStrategy(name, *stratPL, invLib, getRequestHandler());

  // strategy could not be built
  if (strategy == Teuchos::null) {
    RCP<Teuchos::FancyOStream> out = getOutputStream();
    *out << "LSC Construction failed: ";
    *out << "Strategy \"" << name << "\" could not be constructed" << std::endl;
    throw std::runtime_error("LSC Construction failed: Strategy could not be constructed");
  }

  strategy->setSymmetric(isSymmetric_);
  invOpsStrategy_ = strategy;
}

Teuchos::RCP<Teuchos::ParameterList> LSCPreconditionerFactory::getRequestedParameters() const {
  return invOpsStrategy_->getRequestedParameters();
}

bool LSCPreconditionerFactory::updateRequestedParameters(const Teuchos::ParameterList& pl) {
  return invOpsStrategy_->updateRequestedParameters(pl);
}

// Static members and methods

CloneFactory<LSCStrategy> LSCPreconditionerFactory::strategyBuilder_;

RCP<LSCStrategy> LSCPreconditionerFactory::buildStrategy(const std::string& name,
                                                         const Teuchos::ParameterList& settings,
                                                         const RCP<const InverseLibrary>& invLib,
                                                         const RCP<RequestHandler>& rh) {
  Teko_DEBUG_SCOPE("LSCPreconditionerFactory::buildStrategy", 10);

  // initialize the defaults if necessary
  if (strategyBuilder_.cloneCount() == 0) initializeStrategyBuilder();

  // request the preconditioner factory from the CloneFactory
  Teko_DEBUG_MSG("Building LSC strategy \"" << name << "\"", 1);
  RCP<LSCStrategy> strategy = strategyBuilder_.build(name);

  if (strategy == Teuchos::null) return Teuchos::null;

  // now that inverse library has been set,
  // pass in the parameter list
  strategy->setRequestHandler(rh);
  strategy->initializeFromParameterList(settings, *invLib);

  return strategy;
}

void LSCPreconditionerFactory::addStrategy(const std::string& name, const RCP<Cloneable>& clone) {
  // initialize the defaults if necessary
  if (strategyBuilder_.cloneCount() == 0) initializeStrategyBuilder();

  // add clone to builder
  strategyBuilder_.addClone(name, clone);
}

void LSCPreconditionerFactory::initializeStrategyBuilder() {
  RCP<Cloneable> clone;

  // add various strategies to the factory
  clone = rcp(new AutoClone<InvLSCStrategy>());
  strategyBuilder_.addClone("Basic Inverse", clone);

  // add various strategies to the factory
  clone = rcp(new AutoClone<PresLaplaceLSCStrategy>());
  strategyBuilder_.addClone("Pressure Laplace", clone);

  // add various strategies to the factory
  clone = rcp(new AutoClone<LSCSIMPLECStrategy>());
  strategyBuilder_.addClone("SIMPLEC", clone);
}

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