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

bool setDefaultInitialConditionFromNominalValues(
  const Thyra::ModelEvaluator<Scalar>& model,
  const Ptr<StepperBase<Scalar> >& stepper
  )
{

  typedef ScalarTraits<Scalar> ST;
  typedef Thyra::ModelEvaluatorBase MEB;

  if (isInitialized(*stepper))
    return false;  // Already has an initial condition
  
  MEB::InArgs<Scalar> initCond = model.getNominalValues();

  if (!is_null(initCond.get_x())) {
    // IC has x, we will assume that initCont.get_t() is the valid start time.
    // Therefore, we just need to check that x_dot is also set or we will
    // create a zero x_dot
#ifdef RYTHMOS_DEBUG
    THYRA_ASSERT_VEC_SPACES( "setInitialConditionIfExists(...)", 
      *model.get_x_space(), *initCond.get_x()->space() );
#endif
    if (initCond.supports(MEB::IN_ARG_x_dot)) {
      if (is_null(initCond.get_x_dot())) {
        const RCP<Thyra::VectorBase<Scalar> > x_dot =
          createMember(model.get_x_space());
        assign(x_dot.ptr(), ST::zero());
      }
      else {
#ifdef RYTHMOS_DEBUG
        THYRA_ASSERT_VEC_SPACES( "setInitialConditionIfExists(...)", 
          *model.get_x_space(), *initCond.get_x_dot()->space() );
#endif
      }
    }
    stepper->setInitialCondition(initCond);
    return true;
  }

  // The model has not nominal values for which to set the initial
  // conditions so wo don't do anything!  The stepper will still have not
  return false;

}

//.........这里部分代码省略.........
    //
    // This is the linear solve strategy that will be used to solve for the
    // linear system with the W.
    //

    Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
    linearSolverBuilder.setParameterList(sublist(paramList,Stratimikos_name));
    RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >
      W_factory = createLinearSolveStrategy(linearSolverBuilder);

    //
    *out << "\nC) Create and initalize the forward model ...\n";
    //

    // C.1) Create the underlying EpetraExt::ModelEvaluator

    RCP<EpetraExt::DiagonalTransientModel> epetraStateModel =
      EpetraExt::diagonalTransientModel(
        epetra_comm,
        sublist(paramList,DiagonalTransientModel_name)
        );

    *out <<"\nepetraStateModel valid options:\n";
    epetraStateModel->getValidParameters()->print(
      *out, PLPrintOptions().indent(2).showTypes(true).showDoc(true)
      );

    // C.2) Create the Thyra-wrapped ModelEvaluator

    RCP<Thyra::ModelEvaluator<double> > fwdStateModel =
      epetraModelEvaluator(epetraStateModel, W_factory);

    const RCP<const Thyra::VectorSpaceBase<Scalar> >
      x_space = fwdStateModel->get_x_space();

    const RCP<const Thyra::VectorBase<Scalar> >
      gamma = Thyra::create_Vector(epetraStateModel->get_gamma(), x_space);
    *out << "\ngamma = " << describe(*gamma, solnVerbLevel);

    //
    *out << "\nD) Create the stepper and integrator for the forward problem ...\n";
    //

    RCP<Rythmos::TimeStepNonlinearSolver<double> > fwdTimeStepSolver =
      Rythmos::timeStepNonlinearSolver<double>();
    RCP<Rythmos::StepperBase<Scalar> > fwdStateStepper =
      Rythmos::backwardEulerStepper<double>(fwdStateModel, fwdTimeStepSolver);
    fwdStateStepper->setParameterList(sublist(paramList, RythmosStepper_name));
    RCP<Rythmos::IntegratorBase<Scalar> > fwdStateIntegrator;
    {
      RCP<ParameterList>
        integrationControlPL = sublist(paramList, RythmosIntegrationControl_name);
      integrationControlPL->set( "Take Variable Steps", false );
      integrationControlPL->set( "Fixed dt", as<double>(delta_t) );
      RCP<Rythmos::IntegratorBase<Scalar> >
        defaultIntegrator = Rythmos::controlledDefaultIntegrator<Scalar>(
          Rythmos::simpleIntegrationControlStrategy<Scalar>(integrationControlPL)
          );
      fwdStateIntegrator = defaultIntegrator;
    }
    fwdStateIntegrator->setParameterList(sublist(paramList, RythmosIntegrator_name));

    //
    *out << "\nE) Solve the forward problem ...\n";
    //

//.........这里部分代码省略.........
        );
    stateIntegratorAsModel->setVerbLevel(verbLevel);

    *out << "\nUse the StepperAsModelEvaluator to integrate state x(p,finalTime) ... \n";

    RCP<Thyra::VectorBase<Scalar> > x_final;

    {

      Teuchos::OSTab tab(out);

      x_final = createMember(stateIntegratorAsModel->get_g_space(0));

      eval_g(
        *stateIntegratorAsModel,
        0, *state_ic.get_p(0),
        finalTime,
        0, &*x_final
        );

      *out
        << "\nx_final = x(p,finalTime) evaluated using stateIntegratorAsModel:\n"
        << describe(*x_final,solnVerbLevel);

    }

    //
    // Test the integrated state against the exact analytical state solution
    //

    RCP<const Thyra::VectorBase<Scalar> >
      exact_x_final = create_Vector(
        epetraStateModel->getExactSolution(finalTime),
        stateModel->get_x_space()
        );

    result = Thyra::testRelNormDiffErr(
      "exact_x_final", *exact_x_final, "x_final", *x_final,
      "maxStateError", maxStateError, "warningTol", 1.0, // Don't warn
      &*out, solnVerbLevel
      );
    if (!result) success = false;

    //
    // Solve and test the forward sensitivity computation
    //

    if (doFwdSensSolve) {

      //
      // Create the forward sensitivity stepper
      //

      RCP<Rythmos::ForwardSensitivityStepper<Scalar> > stateAndSensStepper =
        Rythmos::forwardSensitivityStepper<Scalar>();
      if (doFwdSensErrorControl) {
        stateAndSensStepper->initializeDecoupledSteppers(
          stateModel, 0, stateModel->getNominalValues(),
          stateStepper, nonlinearSolver,
          integrator->cloneIntegrator(), finalTime
          );
      }
      else {
        stateAndSensStepper->initializeSyncedSteppers(
          stateModel, 0, stateModel->getNominalValues(),
          stateStepper, nonlinearSolver

double computeForwardSensitivityErrorStackedStepperSinCosFE(
    int numTimeSteps,
    Array<RCP<const VectorBase<double> > >& computedSol,
    Array<RCP<const VectorBase<double> > >& exactSol
)
{
    using Teuchos::rcp_dynamic_cast;
    typedef Thyra::ModelEvaluatorBase MEB;
    // Forward ODE Model:
    RCP<SinCosModel> fwdModel = sinCosModel();
    {
        RCP<ParameterList> pl = Teuchos::parameterList();
        pl->set("Accept model parameters",true);
        pl->set("Implicit model formulation",false);
        pl->set("Provide nominal values",true);
        double b = 5.0;
        //double phi = 0.0;
        double a = 2.0;
        double f = 3.0;
        double L = 4.0;
        double x0 = a;
        double x1 = b*f/L;
        pl->set("Coeff a", a);
        pl->set("Coeff f", f);
        pl->set("Coeff L", L);
        pl->set("IC x_0", x0);
        pl->set("IC x_1", x1);
        fwdModel->setParameterList(pl);
    }
    RCP<StepperBase<double> > fwdStepper;
    RCP<StepperBase<double> > fsStepper;
    {
        const RCP<StepperBuilder<double> > builder = stepperBuilder<double>();
        RCP<ParameterList> stepperPL = Teuchos::parameterList();
        stepperPL->set("Stepper Type","Forward Euler");
        builder->setParameterList(stepperPL);
        fwdStepper = builder->create();
        fsStepper = builder->create();
    }
    // Forward Sensitivity Model:
    RCP<ForwardSensitivityExplicitModelEvaluator<double> > fsModel =
        forwardSensitivityExplicitModelEvaluator<double>();
    int p_index = 0;
    fsModel->initializeStructure(fwdModel,p_index);

    const MEB::InArgs<double> fwdModel_ic = fwdModel->getNominalValues();
    fwdStepper->setModel(fwdModel);
    fwdStepper->setInitialCondition(fwdModel_ic);
    fsModel->initializePointState(Teuchos::inOutArg(*fwdStepper),false);

    MEB::InArgs<double> fsModel_ic = fsModel->getNominalValues();
    {
        // Set up sensitivity initial conditions so they match the initial
        // conditions in getExactSensSolution
        RCP<Thyra::VectorBase<double> > s_bar_init
            = createMember(fsModel->get_x_space());
        RCP<Thyra::DefaultMultiVectorProductVector<double> > s_bar_mv =
            rcp_dynamic_cast<Thyra::DefaultMultiVectorProductVector<double> >(
                s_bar_init,
                true
            );
        int np = 3; // SinCos problem number of elements in parameter vector.
        for (int j=0 ; j < np ; ++j) {
            MEB::InArgs<double> sens_ic = fwdModel->getExactSensSolution(j,0.0);
            V_V(outArg(*(s_bar_mv->getNonconstVectorBlock(j))),
                *(sens_ic.get_x())
               );
        }
        fsModel_ic.set_x(s_bar_init);
    }
    fsStepper->setModel(fsModel);
    fsStepper->setInitialCondition(fsModel_ic);

    RCP<StackedStepper<double> > sStepper = stackedStepper<double>();
    sStepper->addStepper(fwdStepper);
    sStepper->addStepper(fsStepper);
    {
        // Set up Forward Sensitivities step strategy
        RCP<ForwardSensitivityStackedStepperStepStrategy<double> > stepStrategy =
            forwardSensitivityStackedStepperStepStrategy<double>();
        sStepper->setStackedStepperStepControlStrategy(stepStrategy);
    }

    double finalTime = 1.0e-4;
    double dt = finalTime/numTimeSteps; // Assume t_0 = 0.0;
    for (int i=0 ; i < numTimeSteps ; ++i ) {
        double dt_taken = sStepper->takeStep(dt,STEP_TYPE_FIXED);
        TEUCHOS_ASSERT( dt_taken == dt );
    }
    RCP<VectorBase<double> > x_bar_final =
        Thyra::createMember(sStepper->get_x_space());
    {
        Array<double> t_vec;
        Array<RCP<const VectorBase<double> > > x_vec;

        t_vec.push_back(finalTime);
        sStepper->getPoints(
            t_vec,
            &x_vec,
            NULL,
//.........这里部分代码省略.........