C++ InArgs::set_x方法代码示例

void
Piro::LOCASolver<Scalar>::evalModelImpl(
    const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs,
    const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs) const
{
  const int l = 0; // TODO: Allow user to select parameter index
  const Teuchos::RCP<const Thyra::VectorBase<Scalar> > p_inargs = inArgs.get_p(l);

  // Forward parameter values to the LOCA stepper
  {
    const Teuchos::RCP<const Thyra::VectorBase<Scalar> > p_inargs_or_nominal =
      Teuchos::nonnull(p_inargs) ? p_inargs : this->getNominalValues().get_p(l);
    const Thyra::ConstDetachedVectorView<Scalar> p_init_values(p_inargs_or_nominal);
    const Teuchos_Ordinal p_entry_count = p_init_values.subDim();
    TEUCHOS_ASSERT(p_entry_count == Teuchos::as<Teuchos_Ordinal>(paramVector_.length()));

    for (Teuchos_Ordinal k = 0; k < p_entry_count; ++k) {
      paramVector_[k] = p_init_values[k];
    }

    group_->setParams(paramVector_);
  }

  stepper_->reset(globalData_, group_, locaStatusTests_, noxStatusTests_, piroParams_);
  const LOCA::Abstract::Iterator::IteratorStatus status = stepper_->run();

  if (status == LOCA::Abstract::Iterator::Finished) {
    std::cerr << "Continuation Stepper Finished.\n";
  } else if (status == LOCA::Abstract::Iterator::NotFinished) {
    std::cerr << "Continuation Stepper did not reach final value.\n";
  } else {
    std::cerr << "Nonlinear solver failed to converge.\n";
    outArgs.setFailed();
  }

  const Teuchos::RCP<Thyra::VectorBase<Scalar> > x_outargs = outArgs.get_g(this->num_g());
  const Teuchos::RCP<Thyra::VectorBase<Scalar> > x_final =
    Teuchos::nonnull(x_outargs) ? x_outargs : Thyra::createMember(this->get_g_space(this->num_g()));

  {
    // Deep copy final solution from LOCA group
    NOX::Thyra::Vector finalSolution(x_final);
    finalSolution = group_->getX();
  }

  // Compute responses for the final solution
  {
    Thyra::ModelEvaluatorBase::InArgs<Scalar> modelInArgs =
      this->getModel().createInArgs();
    {
      modelInArgs.set_x(x_final);
      modelInArgs.set_p(l, p_inargs);
    }

    this->evalConvergedModel(modelInArgs, outArgs);
  }
}

TEUCHOS_UNIT_TEST( Rythmos_ExplicitRKStepper, basePoint ) {
  RCP<SinCosModel> model = sinCosModel(false);
  {
    RCP<ParameterList> pl = Teuchos::parameterList();
    pl->set("Accept model parameters",true);
    model->setParameterList(pl);
  }
  Thyra::ModelEvaluatorBase::InArgs<double> ic = model->getNominalValues();
  // t_ic
  double t_ic = 1.0; // not used
  // x_ic
  RCP<VectorBase<double> > x_ic = Thyra::createMember(*model->get_x_space());
  {
    Thyra::DetachedVectorView<double> x_ic_view( *x_ic );
    x_ic_view[0] = 5.0;
    x_ic_view[1] = 6.0;
  }
  // parameter 0 ic
  RCP<VectorBase<double> > p_ic = Thyra::createMember(*model->get_p_space(0));
  {
    Thyra::DetachedVectorView<double> p_ic_view( *p_ic );
    p_ic_view[0] = 2.0; // a
    p_ic_view[1] = 3.0; // f 
    p_ic_view[2] = 4.0; // L
  }
  ic.set_p(0,p_ic); 
  ic.set_x(x_ic);
  ic.set_t(t_ic);
  RCP<ExplicitRKStepper<double> > stepper = explicitRKStepper<double>();
  stepper->setModel(model);
  stepper->setInitialCondition(ic);
  stepper->setRKButcherTableau(createRKBT<double>("Forward Euler"));
  double dt = 0.2;
  double dt_taken;
  dt_taken = stepper->takeStep(dt,STEP_TYPE_FIXED);
  TEST_EQUALITY_CONST( dt_taken, 0.2 );
  const StepStatus<double> status = stepper->getStepStatus();
  TEST_ASSERT( !is_null(status.solution) );
  double tol = 1.0e-10;
  {
    Thyra::ConstDetachedVectorView<double> x_new_view( *(status.solution) );
    TEST_FLOATING_EQUALITY( x_new_view[0], 5.0 + 0.2*(6.0), tol );
    TEST_FLOATING_EQUALITY( x_new_view[1], 6.0 + 0.2*( (3.0/4.0)*(3.0/4.0)*(2.0-5.0) ), tol );
  }
}

void Piro::RythmosSolver<Scalar>::evalModelImpl(
#endif
    const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs,
    const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs) const
{
    using Teuchos::RCP;
    using Teuchos::rcp;

    // TODO: Support more than 1 parameter and 1 response
    const int j = 0;
    const int l = 0;

    // Parse InArgs
    RCP<const Thyra::VectorBase<Scalar> > p_in;
    if (num_p > 0) {
        p_in = inArgs.get_p(l);
    }
    RCP<const Thyra::VectorBase<Scalar> > p_in2;  //JF add for multipoint
    if (num_p > 1) {
        p_in2 = inArgs.get_p(l+1);
    }

    // Parse OutArgs
    RCP<Thyra::VectorBase<Scalar> > g_out;
    if (num_g > 0) {
        g_out = outArgs.get_g(j);
    }
    const RCP<Thyra::VectorBase<Scalar> > gx_out = outArgs.get_g(num_g);

    Thyra::ModelEvaluatorBase::InArgs<Scalar> state_ic = model->getNominalValues();

    // Set initial time in ME if needed

    if(t_initial > 0.0 && state_ic.supports(Thyra::ModelEvaluatorBase::IN_ARG_t))

        state_ic.set_t(t_initial);

    if (Teuchos::nonnull(initialConditionModel)) {
        // The initial condition depends on the parameter
        // It is found by querying the auxiliary model evaluator as the last response
        const RCP<Thyra::VectorBase<Scalar> > initialState =
            Thyra::createMember(model->get_x_space());

        {
            Thyra::ModelEvaluatorBase::InArgs<Scalar> initCondInArgs = initialConditionModel->createInArgs();
            if (num_p > 0) {
                initCondInArgs.set_p(l, inArgs.get_p(l));
            }

            Thyra::ModelEvaluatorBase::OutArgs<Scalar> initCondOutArgs = initialConditionModel->createOutArgs();
            initCondOutArgs.set_g(initCondOutArgs.Ng() - 1, initialState);

            initialConditionModel->evalModel(initCondInArgs, initCondOutArgs);
        }

        state_ic.set_x(initialState);
    }

    // Set paramters p_in as part of initial conditions
    if (num_p > 0) {
        if (Teuchos::nonnull(p_in)) {
            state_ic.set_p(l, p_in);
        }
    }
    if (num_p > 1) { //JF added for multipoint
        if (Teuchos::nonnull(p_in2)) {
            state_ic.set_p(l+1, p_in2);
        }
    }

    *out << "\nstate_ic:\n" << Teuchos::describe(state_ic, solnVerbLevel);

    //JF  may need a version of the following for multipoint, i.e. num_p>1, l+1, if we want sensitivities
    RCP<Thyra::MultiVectorBase<Scalar> > dgxdp_out;
    Thyra::ModelEvaluatorBase::Derivative<Scalar> dgdp_deriv_out;
    if (num_p > 0) {
        const Thyra::ModelEvaluatorBase::DerivativeSupport dgxdp_support =
            outArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDp, num_g, l);
        if (dgxdp_support.supports(Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM)) {
            const Thyra::ModelEvaluatorBase::Derivative<Scalar> dgxdp_deriv =
                outArgs.get_DgDp(num_g, l);
            dgxdp_out = dgxdp_deriv.getMultiVector();
        }

        if (num_g > 0) {
            const Thyra::ModelEvaluatorBase::DerivativeSupport dgdp_support =
                outArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDp, j, l);
            if (!dgdp_support.none()) {
                dgdp_deriv_out = outArgs.get_DgDp(j, l);
            }
        }
    }

    const bool requestedSensitivities =
        Teuchos::nonnull(dgxdp_out) || !dgdp_deriv_out.isEmpty();

    RCP<const Thyra::VectorBase<Scalar> > finalSolution;
    if (!requestedSensitivities) {
        //
        *out << "\nE) Solve the forward problem ...\n";
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  // The time spent
  globalData_->locaUtils->out() << std::endl <<
    "#### Statistics ########" << std::endl;

  // Check number of steps
  int numSteps = stepper_->getStepNumber();
  globalData_->locaUtils->out() << std::endl <<
    " Number of continuation Steps = " << numSteps << std::endl;

  // Check number of failed steps
  int numFailedSteps = stepper_->getNumFailedSteps();
  globalData_->locaUtils->out() << std::endl <<
    " Number of failed continuation Steps = " << numFailedSteps << std::endl;

  globalData_->locaUtils->out() << std::endl;


  // Note: the last g is used to store the final solution. It can be null - if it is just
  // skip the store. If adaptation has occurred, g is not the correct size.

  const Teuchos::RCP<Thyra::VectorBase<Scalar> > x_outargs = outArgs.get_g(this->num_g());
  Teuchos::RCP<Thyra::VectorBase<Scalar> > x_final;

  int x_args_dim = 0;
  int f_sol_dim = 0;

  // Pardon the nasty cast to resize the last g in outArgs - need to fit the solution
  Thyra::ModelEvaluatorBase::OutArgs<Scalar>* mutable_outArgsPtr =
    const_cast<Thyra::ModelEvaluatorBase::OutArgs<Scalar>* >(&outArgs);

  if(Teuchos::nonnull(x_outargs)){ // g has been allocated, calculate the sizes of g and the solution

    x_args_dim = x_outargs->space()->dim();
//    f_sol_dim = solMgr_->getSolutionGroup()->getX().length();
    f_sol_dim = Teuchos::rcp_dynamic_cast< ::Thyra::LOCAAdaptiveState >(solMgr_->getState())
          ->getSolutionGroup()->getX().length();


  }

  if(Teuchos::is_null(x_outargs) || (x_args_dim != f_sol_dim)){ // g is not the right size

      x_final = Thyra::createMember(this->get_g_space(this->num_g()));

      mutable_outArgsPtr->set_g(this->num_g(), x_final);

  }
  else { // g is OK, use it
    x_final = x_outargs;
  }

  {
    // Deep copy final solution from LOCA group
    NOX::Thyra::Vector finalSolution(x_final);
//    finalSolution = solMgr_->getSolutionGroup()->getX();
    finalSolution = Teuchos::rcp_dynamic_cast< ::Thyra::LOCAAdaptiveState >(solMgr_->getState())
                      ->getSolutionGroup()->getX();

  }

  // If the arrays need resizing
  if(x_args_dim != f_sol_dim){

    const int parameterCount = this->Np();

    for (int pc = 0; pc < parameterCount; ++pc) {
      const Thyra::ModelEvaluatorBase::DerivativeSupport dgdp_support =
        outArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDp, this->num_g(), pc);
      const Thyra::ModelEvaluatorBase::EDerivativeMultiVectorOrientation dgdp_orient =
        Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM;
      if (dgdp_support.supports(dgdp_orient)) {
        const Thyra::ModelEvaluatorBase::DerivativeMultiVector<Scalar> dgdp =
          Thyra::create_DgDp_mv(*this, this->num_g(), pc, dgdp_orient);
        mutable_outArgsPtr->set_DgDp(this->num_g(), pc, dgdp);
      }
    }
  }

  // Compute responses for the final solution
  {
    Thyra::ModelEvaluatorBase::InArgs<Scalar> modelInArgs =
      this->getModel().createInArgs();
    {
      modelInArgs.set_x(x_final);
      modelInArgs.set_p(l, p_inargs);
    }

    this->evalConvergedModel(modelInArgs, outArgs);

    // Save the final solution TODO: this needs to be redone

    Teuchos::RCP<Thyra::ModelEvaluatorBase::InArgs<Scalar> > fp
         = Teuchos::rcp_const_cast<Thyra::ModelEvaluatorBase::InArgs<Scalar> >(finalPoint_);
    Thyra::ModelEvaluatorBase::InArgsSetup<Scalar> ia;
    ia.setSupports(Thyra::ModelEvaluatorBase::IN_ARG_x, true);
    *fp = ia;
    fp->set_x(x_final);

  }
}

//.........这里部分代码省略.........
    double x_1 = dt/(1.0+std::pow(dt*f/L,2))*(-std::pow(f/L,2)*x_ic_0+x_ic_1/dt+std::pow(f/L,2)*a);
    double tol = 1.0e-10;
    Thyra::ConstDetachedVectorView<double> x_final_view( *x_final );
    TEST_FLOATING_EQUALITY( x_final_view[0], x_0, tol );
    TEST_FLOATING_EQUALITY( x_final_view[1], x_1, tol );
  }
  // Copy InterpolationBuffer data into Cubic Spline interpolation buffer for use in Adjoint Solve
  TimeRange<double> fwdTimeRange; 
  RCP<InterpolationBufferBase<double> > fwdCubicSplineInterpBuffer;
  {
    RCP<PointwiseInterpolationBufferAppender<double> > piba = pointwiseInterpolationBufferAppender<double>();
    RCP<InterpolationBuffer<double> > sinkInterpBuffer = interpolationBuffer<double>();
    sinkInterpBuffer->setStorage(storageLimit);
    RCP<CubicSplineInterpolator<double> > csi = cubicSplineInterpolator<double>();
    sinkInterpBuffer->setInterpolator(csi);
    RCP<const InterpolationBufferBase<double> > sourceInterpBuffer;
    {
      RCP<TrailingInterpolationBufferAcceptingIntegratorBase<double> > tibaib = 
        Teuchos::rcp_dynamic_cast<TrailingInterpolationBufferAcceptingIntegratorBase<double> >(fwdIntegrator,true);
      sourceInterpBuffer = tibaib->getTrailingInterpolationBuffer();
    }
    fwdTimeRange = sourceInterpBuffer->getTimeRange();
    piba->append(*sourceInterpBuffer, fwdTimeRange, Teuchos::outArg(*sinkInterpBuffer));
    fwdCubicSplineInterpBuffer = sinkInterpBuffer;

    TimeRange<double> sourceRange = sourceInterpBuffer->getTimeRange();
    TimeRange<double> sinkRange = sinkInterpBuffer->getTimeRange();
    TEST_EQUALITY( sourceRange.lower(), sinkRange.lower() );
    TEST_EQUALITY( sourceRange.upper(), sinkRange.upper() );
  }
  // Adjoint Solve, reading forward solve data from Cubic Spline interpolation buffer
  {
    RCP<ParameterList> ibPL = Teuchos::parameterList();
    ibPL->sublist("Integrator Settings").sublist("Integrator Selection").set("Integrator Type","Default Integrator");
    ibPL->sublist("Integrator Settings").set("Final Time",finalTime);
    ibPL->sublist("Integration Control Strategy Selection").set("Integration Control Strategy Type","Simple Integration Control Strategy");
    ibPL->sublist("Integration Control Strategy Selection").sublist("Simple Integration Control Strategy").set("Take Variable Steps",false);
    ibPL->sublist("Integration Control Strategy Selection").sublist("Simple Integration Control Strategy").set("Fixed dt",dt);

    ibPL->sublist("Stepper Settings").sublist("Stepper Selection").set("Stepper Type","Backward Euler");
    //ibPL->sublist("Stepper Settings").sublist("Stepper Selection").set("Stepper Type","Implicit RK");
    //ibPL->sublist("Stepper Settings").sublist("Runge Kutta Butcher Tableau Selection").set("Runge Kutta Butcher Tableau Type","Implicit 1 Stage 2nd order Gauss");
    ibPL->sublist("Interpolation Buffer Settings").sublist("Trailing Interpolation Buffer Selection").set("Interpolation Buffer Type","Interpolation Buffer");
    ibPL->sublist("Interpolation Buffer Settings").sublist("Trailing Interpolation Buffer Selection").sublist("Interpolation Buffer").set("StorageLimit",storageLimit);
    ibPL->sublist("Interpolation Buffer Settings").sublist("Interpolator Selection").set("Interpolator Type","Linear Interpolator");
    ib->setParameterList(ibPL);
  }
  RCP<Thyra::ModelEvaluator<double> > adjModel;
  {
    RCP<Rythmos::AdjointModelEvaluator<double> > model = 
      Rythmos::adjointModelEvaluator<double>(
          fwdModel, fwdTimeRange
          );
    //model->setFwdStateSolutionBuffer(fwdCubicSplineInterpBuffer);
    adjModel = model;
  }
  Thyra::ModelEvaluatorBase::InArgs<double> adjIC = adjModel->getNominalValues();
  double phi_ic_0 = 2.0;
  double phi_ic_1 = 3.0;
  {
    // Initial conditions for adjoint:
    const RCP<const Thyra::VectorSpaceBase<double> >
      f_space = fwdModel->get_f_space();
    const RCP<Thyra::VectorBase<double> > x_ic = createMember(f_space);
    {
      Thyra::DetachedVectorView<double> x_ic_view( *x_ic );
      x_ic_view[0] = phi_ic_0;
      x_ic_view[1] = phi_ic_1;
    }
    const RCP<Thyra::VectorBase<double> > xdot_ic = createMember(f_space);
    V_S( Teuchos::outArg(*xdot_ic), ST::zero() );
    adjIC.set_x(x_ic);
    adjIC.set_x_dot(xdot_ic);
  }
  RCP<Thyra::LinearNonlinearSolver<double> > adjNLSolver = Thyra::linearNonlinearSolver<double>();
  RCP<IntegratorBase<double> > adjIntegrator = ib->create(adjModel,adjIC,adjNLSolver);
  RCP<const VectorBase<double> > phi_final;
  {
    Array<double> time_vec;
    time_vec.push_back(finalTime);
    Array<RCP<const Thyra::VectorBase<double> > > phi_final_array;
    adjIntegrator->getFwdPoints(time_vec,&phi_final_array,NULL,NULL);
    phi_final = phi_final_array[0];
  }
  // Verify phi_final is correct
  {
    // Defaults from SinCos Model:
    double f = 1.0;
    double L = 1.0;
    double h = -dt;
    double phi_0 = 1.0/(1.0+std::pow(f*h/L,2.0))*(phi_ic_0+std::pow(f/L,2.0)*h*phi_ic_1);
    double phi_1 = 1.0/(1.0+std::pow(f*h/L,2.0))*(-h*phi_ic_0+phi_ic_1);
    double tol = 1.0e-10;
    Thyra::ConstDetachedVectorView<double> phi_final_view( *phi_final );
    TEST_FLOATING_EQUALITY( phi_final_view[0], phi_0, tol );
    TEST_FLOATING_EQUALITY( phi_final_view[1], phi_1, tol );
  }
  // Compute error estimate
  //TEST_ASSERT( false );
}

void
Piro::MatrixFreeLinearOp<Scalar>::applyImpl(
    const Thyra::EOpTransp M_trans,
    const Thyra::MultiVectorBase<Scalar> &X,
    const Teuchos::Ptr<Thyra::MultiVectorBase<Scalar> > &Y,
    const Scalar alpha,
    const Scalar beta) const
{
  using Teuchos::RCP;
  using Teuchos::Ptr;

  TEUCHOS_TEST_FOR_EXCEPTION(
      !this->opSupported(M_trans),
      Thyra::Exceptions::OpNotSupported,
      this->description() << " does not support operation " << Thyra::toString(M_trans));

  TEUCHOS_TEST_FOR_EXCEPTION(
      !X.range()->isCompatible(*this->domain()),
      Thyra::Exceptions::IncompatibleVectorSpaces,
      "Domain of " << this->description() << ": " << this->domain()->description() <<
      " is not compatible with column space of " << X.description() << ": " << X.range()->description());

  TEUCHOS_TEST_FOR_EXCEPTION(
      !Y->range()->isCompatible(*this->range()),
      Thyra::Exceptions::IncompatibleVectorSpaces,
      "Range of " << this->description() << ": " << this->range()->description() <<
      " is not compatible with column space of " << Y->description() << ": " << Y->range()->description());

  TEUCHOS_TEST_FOR_EXCEPTION(
      !Y->domain()->isCompatible(*X.domain()),
      Thyra::Exceptions::IncompatibleVectorSpaces,
      "Row space of " << Y->description() << ": " << Y->domain()->description() <<
      " is not compatible with row space of " << X.description() << ": " << X.domain()->description());

  TEUCHOS_TEST_FOR_EXCEPTION(
      &X == Y.get(),
      std::logic_error,
      "X and Y arguments are both aliases of " << X.description());

  if (alpha == Teuchos::ScalarTraits<Scalar>::zero()) {
    // Y <- beta * Y
    Thyra::Vt_S(Y, beta);
    return;
  }

  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMagnitude;

  RCP<const Thyra::VectorBase<Scalar> > x_dot_base;
  if (basePoint_.supports(Thyra::ModelEvaluatorBase::IN_ARG_x_dot)) 
    x_dot_base = basePoint_.get_x_dot();

  RCP<const Thyra::VectorBase<Scalar> > x_base = basePoint_.get_x();
  if (Teuchos::is_null(x_base)) {
    x_base = model_->getNominalValues().get_x();
  }
  x_base.assert_not_null();

  const ScalarMagnitude norm_x_base = Thyra::norm_2(*x_base);

  // Number of columns common to both vectors X and Y
  // (X and Y have compatible row spaces)
  const Thyra::Ordinal colCount = X.domain()->dim();
  for (Teuchos::Ordinal j = Teuchos::Ordinal(); j < colCount; ++j) {
    const RCP<const Thyra::VectorBase<Scalar> > X_vec = X.col(j);
    const RCP<Thyra::VectorBase<Scalar> > Y_vec = Y->col(j);

    const ScalarMagnitude norm_dx = Thyra::norm_2(*X_vec);

    if (norm_dx == Teuchos::ScalarTraits<ScalarMagnitude>::zero()) {
      if (beta == Teuchos::ScalarTraits<Scalar>::zero()) {
        // Y_vec <- 0
        Thyra::put_scalar(Teuchos::ScalarTraits<ScalarMagnitude>::zero(), Y_vec.ptr());
      } else {
        // Y_vec <- beta * Y_vec
        Thyra::scale(beta, Y_vec.ptr());
      }
    } else {
      // Scalar perturbation
      const ScalarMagnitude relative_pert_ratio = static_cast<ScalarMagnitude>(lambda_);
      const ScalarMagnitude eta = (relative_pert_ratio * ((norm_x_base / norm_dx) + relative_pert_ratio));

      // Compute perturbed residual
      // Dynamic: f_pert <- f(x_dot_base + eta * (W_alpha * X), x_base + eta * (W_beta * X))
      // Static: f_pert <- f(x_base + eta * X)
      const RCP<Thyra::VectorBase<Scalar> > f_pert = Thyra::createMember(this->range());
      {
        Thyra::ModelEvaluatorBase::InArgs<Scalar> pertInArgs = model_->createInArgs();
        {
          pertInArgs.setArgs(basePoint_);

          const bool isDynamic = Teuchos::nonnull(x_dot_base);

          if (isDynamic) {
            const RCP<Thyra::VectorBase<Scalar> > x_dot_pert = Thyra::createMember(this->domain());
            const Scalar W_alpha = pertInArgs.get_alpha();
            Thyra::V_VpStV<Scalar>(x_dot_pert.ptr(), *x_dot_base, W_alpha * eta, *X_vec);
            pertInArgs.set_x_dot(x_dot_pert);
          }

          const RCP<Thyra::VectorBase<Scalar> > x_pert = Thyra::createMember(this->domain());
//.........这里部分代码省略.........