Tempus_SolutionState_impl.hpp

Go to the documentation of this file.

// @HEADER
// ****************************************************************************
//                Tempus: Copyright (2017) Sandia Corporation
//
// Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
// ****************************************************************************
// @HEADER

#ifndef Tempus_SolutionState_impl_hpp
#define Tempus_SolutionState_impl_hpp

#include "Thyra_VectorStdOps.hpp"

namespace Tempus {


template<class Scalar>
SolutionState<Scalar>::SolutionState()
  : x_              (Teuchos::null),
    x_nc_           (Teuchos::null),
    xdot_           (Teuchos::null),
    xdot_nc_        (Teuchos::null),
    xdotdot_        (Teuchos::null),
    xdotdot_nc_     (Teuchos::null),
    stepperState_   (Teuchos::null),
    stepperState_nc_(Teuchos::null),
    physicsState_   (Teuchos::null),
    physicsState_nc_(Teuchos::null)
{
  metaData_nc_     = Teuchos::rcp(new SolutionStateMetaData<Scalar>());
  metaData_        = metaData_nc_;
  stepperState_nc_ = Teuchos::rcp(new StepperState<Scalar>("Default"));
  stepperState_    = stepperState_nc_;
  physicsState_nc_ = Teuchos::rcp(new PhysicsState<Scalar> ());
  physicsState_    = physicsState_nc_;
}


template<class Scalar>
SolutionState<Scalar>::SolutionState(
  const Teuchos::RCP<SolutionStateMetaData<Scalar> > metaData,
  const Teuchos::RCP<Thyra::VectorBase<Scalar> >& x,
  const Teuchos::RCP<Thyra::VectorBase<Scalar> >& xdot,
  const Teuchos::RCP<Thyra::VectorBase<Scalar> >& xdotdot,
  const Teuchos::RCP<StepperState<Scalar> >& stepperState,
  const Teuchos::RCP<PhysicsState<Scalar> >& physicsState)
  : metaData_       (metaData),
    metaData_nc_    (metaData),
    x_              (x),
    x_nc_           (x),
    xdot_           (xdot),
    xdot_nc_        (xdot),
    xdotdot_        (xdotdot),
    xdotdot_nc_     (xdotdot),
    stepperState_   (stepperState),
    stepperState_nc_(stepperState),
    physicsState_   (physicsState),
    physicsState_nc_(physicsState)
{
  if (stepperState_nc_ == Teuchos::null) {
    stepperState_nc_ = Teuchos::rcp(new StepperState<Scalar>("Default"));
    stepperState_    = stepperState_nc_;
  }
  if (physicsState_nc_ == Teuchos::null) {
    physicsState_nc_ = Teuchos::rcp(new PhysicsState<Scalar> ());
    physicsState_    = physicsState_nc_;
  }
}

template<class Scalar>
SolutionState<Scalar>::SolutionState(
  const Teuchos::RCP<const SolutionStateMetaData<Scalar> > metaData,
  const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& x,
  const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& xdot,
  const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& xdotdot,
  const Teuchos::RCP<const StepperState<Scalar> >& stepperState,
  const Teuchos::RCP<const PhysicsState<Scalar> >& physicsState)
  : metaData_       (metaData),
    metaData_nc_    (Teuchos::null),
    x_              (x),
    x_nc_           (Teuchos::null),
    xdot_           (xdot),
    xdot_nc_        (Teuchos::null),
    xdotdot_        (xdotdot),
    xdotdot_nc_     (Teuchos::null),
    stepperState_   (stepperState),
    stepperState_nc_(Teuchos::null),
    physicsState_   (physicsState),
    physicsState_nc_(Teuchos::null)
{
  if (stepperState_ == Teuchos::null) {
    stepperState_ = Teuchos::rcp(new StepperState<Scalar>("Default"));
  }
  if (physicsState_ == Teuchos::null) {
    physicsState_ = Teuchos::rcp(new PhysicsState<Scalar> ());
  }
}


template<class Scalar>
SolutionState<Scalar>::SolutionState(const SolutionState<Scalar>& ss_)
  :metaData_       (ss_.metaData_),
   metaData_nc_    (ss_.metaData_nc_),
   x_              (ss_.x_),
   x_nc_           (ss_.x_nc_),
   xdot_           (ss_.xdot_),
   xdot_nc_        (ss_.xdot_nc_),
   xdotdot_        (ss_.xdotdot_),
   xdotdot_nc_     (ss_.xdotdot_nc_),
   stepperState_   (ss_.stepperState_),
   stepperState_nc_(ss_.stepperState_nc_),
   physicsState_   (ss_.physicsState_),
   physicsState_nc_(ss_.physicsState_nc_)
{}


template<class Scalar>
Teuchos::RCP<SolutionState<Scalar> > SolutionState<Scalar>::clone() const
{
  using Teuchos::RCP;

  RCP<SolutionStateMetaData<Scalar> > metaData_out;
  if (!Teuchos::is_null(metaData_)) metaData_out = metaData_->clone();

  RCP<Thyra::VectorBase<Scalar> > x_out;
  if (!Teuchos::is_null(x_)) x_out = x_->clone_v();

  RCP<Thyra::VectorBase<Scalar> > xdot_out;
  if (!Teuchos::is_null(xdot_)) xdot_out = xdot_->clone_v();

  RCP<Thyra::VectorBase<Scalar> > xdotdot_out;
  if (!Teuchos::is_null(xdotdot_)) xdotdot_out = xdotdot_->clone_v();

  RCP<StepperState<Scalar> > sS_out;
  if (!Teuchos::is_null(stepperState_)) sS_out=stepperState_->clone();

  RCP<PhysicsState<Scalar> > pS_out;
  if (!Teuchos::is_null(physicsState_)) pS_out=physicsState_->clone();

  RCP<SolutionState<Scalar> > ss_out = Teuchos::rcp(new SolutionState<Scalar> (
    metaData_out, x_out, xdot_out, xdotdot_out, sS_out, pS_out));

  return ss_out;
}


template<class Scalar>
void SolutionState<Scalar>::
copy(const Teuchos::RCP<const SolutionState<Scalar> >& ss)
{
  metaData_nc_->copy(ss->metaData_);
  this->copySolutionData(ss);
}


template<class Scalar>
void SolutionState<Scalar>::
copySolutionData(const Teuchos::RCP<const SolutionState<Scalar> >& ss)
{
  if (ss->x_ == Teuchos::null)
    x_nc_ = Teuchos::null;
  else {
    if (x_nc_ == Teuchos::null) {
      x_nc_ = ss->x_->clone_v();
    }
    else
      Thyra::V_V(x_nc_.ptr(),       *(ss->x_));
  }
  x_ = x_nc_;

  if (ss->xdot_ == Teuchos::null)
    xdot_nc_ = Teuchos::null;
  else {
    if (xdot_nc_ == Teuchos::null)
      xdot_nc_ = ss->xdot_->clone_v();
    else
      Thyra::V_V(xdot_nc_.ptr(),    *(ss->xdot_));
  }
  xdot_ = xdot_nc_;

  if (ss->xdotdot_ == Teuchos::null)
    xdotdot_nc_ = Teuchos::null;
  else {
    if (xdotdot_nc_ == Teuchos::null)
      xdotdot_nc_ = ss->xdotdot_->clone_v();
    else
      Thyra::V_V(xdotdot_nc_.ptr(), *(ss->xdotdot_));
  }
  xdotdot_ = xdotdot_nc_;

  if (ss->stepperState_ == Teuchos::null)
    stepperState_nc_ = Teuchos::null;
  else {
    if (stepperState_nc_ == Teuchos::null)
      stepperState_nc_ = ss->stepperState_->clone();
    else
      stepperState_nc_->copy(ss->stepperState_);
  }
  stepperState_ = stepperState_nc_;

  if (ss->physicsState_ == Teuchos::null)
    physicsState_nc_ = Teuchos::null;
  else {
    if (physicsState_nc_ == Teuchos::null)
      physicsState_nc_ = ss->physicsState_->clone();
    else
      physicsState_nc_->copy(ss->physicsState_);
  }
  physicsState_ = physicsState_nc_;
}

template<class Scalar>
bool SolutionState<Scalar>::operator< (const SolutionState<Scalar>& ss) const
{
  return (this->metaData_->getTime() < ss.metaData_->getTime());
}

template<class Scalar>
bool SolutionState<Scalar>::operator<= (const SolutionState<Scalar>& ss) const
{
  return (this->metaData_->getTime() <= ss.metaData_->getTime());
}

template<class Scalar>
bool SolutionState<Scalar>::operator< (const Scalar& t) const
{
  return (this->metaData_->getTime() < t);
}

template<class Scalar>
bool SolutionState<Scalar>::operator<= (const Scalar& t) const
{
  return (this->metaData_->getTime() <= t);
}

template<class Scalar>
bool SolutionState<Scalar>::operator> (const SolutionState<Scalar>& ss) const
{
  return (this->metaData_->getTime() > ss.metaData_->getTime());
}

template<class Scalar>
bool SolutionState<Scalar>::operator>= (const SolutionState<Scalar>& ss) const
{
  return (this->metaData_->getTime() >= ss.metaData_->getTime());
}

template<class Scalar>
bool SolutionState<Scalar>::operator> (const Scalar& t) const
{
  return (this->metaData_->getTime() > t);
}

template<class Scalar>
bool SolutionState<Scalar>::operator>= (const Scalar& t) const
{
  return (this->metaData_->getTime() >= t);
}

template<class Scalar>
bool SolutionState<Scalar>::operator== (const SolutionState<Scalar>& ss) const
{
  return (this->metaData_->getTime() == ss.metaData_->getTime());
}

template<class Scalar>
bool SolutionState<Scalar>::operator== (const Scalar& t) const
{
  return (this->metaData_->getTime() == t);
}

template<class Scalar>
std::string SolutionState<Scalar>::description() const
{
  std::ostringstream out;
  out << "SolutionState"
      << " (index =" <<std::setw(6)<< this->getIndex()
      << "; time =" <<std::setw(10)<<std::setprecision(3)<<this->getTime()
      << "; dt ="   <<std::setw(10)<<std::setprecision(3)<<this->getTimeStep()
      << ")";
  return out.str();
}

template<class Scalar>
void SolutionState<Scalar>::describe(
   Teuchos::FancyOStream               &out,
   const Teuchos::EVerbosityLevel      verbLevel) const
{
  auto l_out = Teuchos::fancyOStream( out.getOStream() );
  Teuchos::OSTab ostab(*l_out, 2, this->description());
  l_out->setOutputToRootOnly(0);

  *l_out << "\n--- " << this->description() << " ---" << std::endl;

  if (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_EXTREME)) {

    metaData_->describe(*l_out,verbLevel);
    *l_out << "  x       = " << std::endl;
    x_->describe(*l_out,verbLevel);

    if (xdot_ != Teuchos::null) {
      *l_out << "  xdot_   = " << std::endl;
      xdot_->describe(*l_out,verbLevel);
    }
    if (xdotdot_ != Teuchos::null) {
      *l_out << "  xdotdot = " << std::endl;
      xdotdot_->describe(*l_out,verbLevel);
    }

    if (stepperState_ != Teuchos::null)
      stepperState_->describe(*l_out,verbLevel);
    if (physicsState_ != Teuchos::null)
      physicsState_->describe(*l_out,verbLevel);

    *l_out << std::string(this->description().length()+8, '-') <<std::endl;
  }
}


template<class Scalar>
void SolutionState<Scalar>::computeNorms(
  const Teuchos::RCP<const SolutionState<Scalar> >& ssIn)
{
  if (!getComputeNorms()) return;

  auto x   = this->getX();
  this->setXNormL2(Thyra::norm(*x));

  if (ssIn != Teuchos::null) {
    auto xIn = ssIn->getX();

    // dx = x - xIn
    Teuchos::RCP<Thyra::VectorBase<Scalar> > dx = Thyra::createMember(x->space());
    Thyra::V_VmV(dx.ptr(), *x, *xIn);
    Scalar dxNorm = Thyra::norm(*dx);
    Scalar xInNorm = Thyra::norm(*xIn);
    this->setDxNormL2Abs(dxNorm);
    // Compute change, e.g., ||x^n-x^(n-1)||/||x^(n-1)||
    const Scalar eps = std::numeric_limits<Scalar>::epsilon();
    const Scalar min = std::numeric_limits<Scalar>::min();
    if ( xInNorm < min/eps ) {  // numerically zero
      this->setDxNormL2Rel(std::numeric_limits<Scalar>::infinity());
    } else {
      //this->setDxNormL2Rel(dxNorm/(xInNorm + eps));
      this->setDxNormL2Rel(dxNorm/(xInNorm*(1.0 + 1.0e4*eps)));
    }
  }
}


// Nonmember constructors.
// ------------------------------------------------------------------------

template<class Scalar>
Teuchos::RCP<SolutionState<Scalar> > createSolutionStateX(
  const Teuchos::RCP<Thyra::VectorBase<Scalar> >& x,
  const Teuchos::RCP<Thyra::VectorBase<Scalar> >& xdot,
  const Teuchos::RCP<Thyra::VectorBase<Scalar> >& xdotdot)
{
  Teuchos::RCP<SolutionStateMetaData<Scalar> >
    metaData_nc = Teuchos::rcp(new SolutionStateMetaData<Scalar>());

  Teuchos::RCP<StepperState<Scalar> >
    stepperState_nc = Teuchos::rcp(new StepperState<Scalar>("Default"));

  Teuchos::RCP<PhysicsState<Scalar> >
    physicsState_nc = Teuchos::rcp(new PhysicsState<Scalar> ());

  Teuchos::RCP<SolutionState<Scalar> > ss = rcp(new SolutionState<Scalar> (
    metaData_nc, x, xdot, xdotdot, stepperState_nc, physicsState_nc));

  return ss;
}

template<class Scalar>
Teuchos::RCP<SolutionState<Scalar> > createSolutionStateX(
  const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& x,
  const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& xdot,
  const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& xdotdot)
{
  Teuchos::RCP<const SolutionStateMetaData<Scalar> >
    metaData = Teuchos::rcp(new SolutionStateMetaData<Scalar>());

  Teuchos::RCP<const StepperState<Scalar> >
    stepperState = Teuchos::rcp(new StepperState<Scalar>("Default"));

  Teuchos::RCP<const PhysicsState<Scalar> >
    physicsState = Teuchos::rcp(new PhysicsState<Scalar> ());

  Teuchos::RCP<SolutionState<Scalar> > ss = rcp(new SolutionState<Scalar> (
    metaData, x, xdot, xdotdot, stepperState, physicsState));

  return ss;
}


template<class Scalar>
Teuchos::RCP<SolutionState<Scalar> > createSolutionStateME(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
  const Teuchos::RCP<StepperState<Scalar> >& stepperState,
  const Teuchos::RCP<PhysicsState<Scalar> >& physicsState)
{
  typedef Thyra::ModelEvaluatorBase MEB;
  using Teuchos::rcp_const_cast;

  auto metaData_nc = Teuchos::rcp(new SolutionStateMetaData<Scalar>());
  metaData_nc->setSolutionStatus(Status::PASSED);

  MEB::InArgs<Scalar> inArgs = model->getNominalValues();

  TEUCHOS_TEST_FOR_EXCEPTION(inArgs.supports(MEB::IN_ARG_x) == false,
    std::logic_error,
    model->description() << "does not support an x solution vector!");

  // The solution vector, x, is required (usually).
  auto x_nc = rcp_const_cast<Thyra::VectorBase<Scalar> > (inArgs.get_x());

  // The solution derivative, xdot, can be optional provided, based on
  // application needs.  Here we will base it on "supports" IN_ARG_x_dot.
  // Depending on the stepper used, a temporary xdot vector may be created
  // within the Stepper, but not moved to the SolutionState.
  Teuchos::RCP<Thyra::VectorBase<Scalar> > xdot_nc;
  if (inArgs.supports(MEB::IN_ARG_x_dot)) {
    xdot_nc = rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs.get_x_dot());
  } else {
    xdot_nc = Teuchos::null;
  }

  // Similar as xdot.
  Teuchos::RCP<Thyra::VectorBase<Scalar> > xdotdot_nc;
  if (inArgs.supports(MEB::IN_ARG_x_dot_dot)) {
    xdotdot_nc =
      rcp_const_cast<Thyra::VectorBase<Scalar> > (inArgs.get_x_dot_dot());
  } else {
    xdotdot_nc = Teuchos::null;
  }

  Teuchos::RCP<StepperState<Scalar> > stepperState_nc;
  if (stepperState == Teuchos::null) {
    stepperState_nc = Teuchos::rcp(new StepperState<Scalar> ()); // Use default
  } else {
    stepperState_nc = stepperState;
  }

  Teuchos::RCP<PhysicsState<Scalar> > physicsState_nc;
  if (physicsState == Teuchos::null) {
    physicsState_nc = Teuchos::rcp(new PhysicsState<Scalar> ()); // Use default
  } else {
    physicsState_nc = physicsState;
  }

  Teuchos::RCP<SolutionState<Scalar> > ss = rcp(new SolutionState<Scalar> (
    metaData_nc, x_nc, xdot_nc, xdotdot_nc, stepperState_nc, physicsState_nc));

  return ss;
}

} // namespace Tempus
#endif // Tempus_SolutionState_impl_hpp