C++ TimeRange::upper方法代码示例

TEUCHOS_UNIT_TEST( Rythmos_TimeRange, copyAndScaleInvalid ) {
  TimeRange<double> tr;
  TimeRange<double> newTr = tr.copyAndScale(5.0);
  TEST_EQUALITY_CONST( newTr.isValid(), false );
  TEST_EQUALITY( newTr.lower(), tr.lower() );
  TEST_EQUALITY( newTr.upper(), tr.upper() );
  TEST_EQUALITY( newTr.length(), tr.length() );
}

bool Rythmos::isInRange_cc(const TimeRange<TimeType> &tr, const TimeType &p)
{
  return (
    compareTimeValues(p,tr.lower()) >= 0
    && compareTimeValues(p,tr.upper()) <= 0
    );
}

TEUCHOS_UNIT_TEST( Rythmos_TimeRange, copyAndScale ) {
  TimeRange<double> tr(1.0,2.0);
  TimeRange<double> newTr = tr.copyAndScale(5.0);
  TEST_EQUALITY_CONST( newTr.isValid(), true );
  TEST_EQUALITY_CONST( newTr.lower(), 5.0 );
  TEST_EQUALITY_CONST( newTr.upper(), 10.0 );
  TEST_EQUALITY_CONST( newTr.length(), 5.0 );
}

void Rythmos::assertNoTimePointsInsideCurrentTimeRange(
  const InterpolationBufferBase<Scalar>& interpBuffer,
  const Array<Scalar>& time_vec
  )
{
  typedef ScalarTraits<Scalar> ST;
  const int numTimePoints = time_vec.size();
  const TimeRange<Scalar> currentTimeRange = interpBuffer.getTimeRange();
  if (currentTimeRange.length() >= ST::zero()) {
    for ( int i = 0; i < numTimePoints; ++i ) {
      TEST_FOR_EXCEPTION(
        currentTimeRange.isInRange(time_vec[i]), std::out_of_range,
        "Error, time_vec["<<i<<"] = " << time_vec[i] << " is in TimeRange of " 
        << interpBuffer.description() << " = ["
        << currentTimeRange.lower() << "," << currentTimeRange.upper() << "]!"
        );
    }
  }
}

TEUCHOS_UNIT_TEST( BasicDiscreteAdjointStepperTester, rawNonlinearAdjoint )
{

  using Teuchos::outArg;
  using Teuchos::describe;
  using Teuchos::getParametersFromXmlString;
  typedef Thyra::ModelEvaluatorBase MEB;

  //
  out << "\nA) Create the nonlinear ME ...\n";
  //

  RCP<VanderPolModel> stateModel = vanderPolModel(
    getParametersFromXmlString(
      "<ParameterList>"
      "  <Parameter name=\"Implicit model formulation\" type=\"bool\" value=\"1\"/>"
      "</ParameterList>"
      )
    );

  //
  out << "\nB) Create the nonlinear solver ...\n";
  //

  RCP<TimeStepNonlinearSolver<double> > nlSolver = timeStepNonlinearSolver<double>(
    getParametersFromXmlString(
      "<ParameterList>"
      "  <Parameter name=\"Default Tol\" type=\"double\" value=\"1.0e-10\"/>"
      "  <Parameter name=\"Default Max Iters\" type=\"int\" value=\"20\"/>"
      "</ParameterList>"
      )
    );

  //
  out << "\nC) Create the integrator for the forward state problem ...\n";
  //

  RCP<IntegratorBuilder<double> > ib = integratorBuilder<double>(
    Teuchos::getParametersFromXmlString(
      "<ParameterList>"
      "  <ParameterList name=\"Stepper Settings\">"
      "    <ParameterList name=\"Stepper Selection\">"
      "      <Parameter name=\"Stepper Type\" type=\"string\" value=\"Backward Euler\"/>"
      "    </ParameterList>"
      "  </ParameterList>"
      "  <ParameterList name=\"Integration Control Strategy Selection\">"
      "    <Parameter name=\"Integration Control Strategy Type\" type=\"string\""
      "      value=\"Simple Integration Control Strategy\"/>"
      "    <ParameterList name=\"Simple Integration Control Strategy\">"
      "      <Parameter name=\"Take Variable Steps\" type=\"bool\" value=\"false\"/>"
      "      <Parameter name=\"Fixed dt\" type=\"double\" value=\"0.5\"/>" // Gives 2 time steps!
      "    </ParameterList>"
      "  </ParameterList>"
      "  <ParameterList name=\"Interpolation Buffer Settings\">"
      "    <ParameterList name=\"Trailing Interpolation Buffer Selection\">"
      "      <Parameter name=\"Interpolation Buffer Type\" type=\"string\" value=\"Interpolation Buffer\"/>"
      "    </ParameterList>"
      "  </ParameterList>"
      "</ParameterList>"
      )
    );
  
  MEB::InArgs<double> ic = stateModel->getNominalValues();
  RCP<IntegratorBase<double> > integrator = ib->create(stateModel, ic, nlSolver);
  //integrator->setVerbLevel(Teuchos::VERB_EXTREME);

  // ToDo: Set the trailing IB to pick up the entire state solution!

  // 
  out << "\nD) Solve the basic forward problem ...\n";
  //

  const TimeRange<double> fwdTimeRange = integrator->getFwdTimeRange();
  const double t_final = fwdTimeRange.upper();
  RCP<const Thyra::VectorBase<double> > x_final, x_dot_final;
  get_fwd_x_and_x_dot( *integrator, t_final, outArg(x_final), outArg(x_dot_final) );

  out << "\nt_final = " << t_final << "\n";
  out << "\nx_final: " << *x_final;
  out << "\nx_dot_final: " << *x_dot_final;

  //
  out << "\nE) Create the basic adjoint model (no distributed response) ...\n";
  //

  RCP<AdjointModelEvaluator<double> > adjModel =
    adjointModelEvaluator<double>(
      stateModel, fwdTimeRange
      );
  adjModel->setFwdStateSolutionBuffer(integrator);

  //
  out << "\nF) Create a stepper and integrator for the adjoint ...\n";
  //
  
  RCP<Thyra::LinearNonlinearSolver<double> > adjTimeStepSolver =
    Thyra::linearNonlinearSolver<double>();
  RCP<Rythmos::StepperBase<double> > adjStepper =
    integrator->getStepper()->cloneStepperAlgorithm();

//.........这里部分代码省略.........

void PointwiseInterpolationBufferAppender<Scalar>::append(
  const InterpolationBufferBase<Scalar>& interpBuffSource, 
  const TimeRange<Scalar>& appendRange,
  const Ptr<InterpolationBufferBase<Scalar> > &interpBuffSink 
  ) 
{
  TEUCHOS_ASSERT( !is_null(interpBuffSink) );
#ifdef RYTHMOS_DEBUG
  this->assertAppendPreconditions(interpBuffSource,appendRange,*interpBuffSink);
#endif // RYTHMOS_DEBUG

  RCP<Teuchos::FancyOStream> out = this->getOStream();
  Teuchos::OSTab ostab(out,1,"PointwiseInterpolationBufferAppender::append");
  if ( Teuchos::as<int>(this->getVerbLevel()) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
    *out << "Interpolation Buffer source range = [" << interpBuffSource.getTimeRange().lower() << "," <<
      interpBuffSource.getTimeRange().upper() << "]" << std::endl;
    *out << "Append range = [" << appendRange.lower() << "," << appendRange.upper() << "]" << std::endl;
    *out << "Interpolation Buffer sink range = [" << interpBuffSink->getTimeRange().lower() << "," <<
      interpBuffSink->getTimeRange().upper() << "]" << std::endl;
  }
  // Set up appendRange correctly to be either (] or [):
  RCP<const TimeRange<Scalar> > correctedAppendRange = Teuchos::rcp(&appendRange,false);
  if (compareTimeValues<Scalar>(interpBuffSink->getTimeRange().upper(),appendRange.lower()) == 0) {
    // adding to end of buffer 
    correctedAppendRange = Teuchos::rcp(new TimeRange_oc<Scalar>(appendRange));
    if ( Teuchos::as<int>(this->getVerbLevel()) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
      *out << "Corrected append range = (" << correctedAppendRange->lower() << "," << 
        correctedAppendRange->upper() << "]" << std::endl;
    }
  } 
  else if (compareTimeValues<Scalar>(interpBuffSink->getTimeRange().lower(),appendRange.upper()) == 0) {
    // adding to beginning of buffer
    correctedAppendRange = Teuchos::rcp(new TimeRange_co<Scalar>(appendRange));
    if ( Teuchos::as<int>(this->getVerbLevel()) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
      *out << "Corrected append range = [" << correctedAppendRange->lower() << "," << 
        correctedAppendRange->upper() << ")" << std::endl;
    }
  }

  Array<Scalar> time_vec_in;
  interpBuffSource.getNodes(&time_vec_in);

  Array<Scalar> time_vec;
  selectPointsInTimeRange(time_vec_in,*correctedAppendRange,Teuchos::outArg(time_vec));
  if ( Teuchos::as<int>(this->getVerbLevel()) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
    *out << "Selected points for appending to sink buffer: " << time_vec << std::endl;
  }

  Array<RCP<const Thyra::VectorBase<Scalar> > > x_vec;
  Array<RCP<const Thyra::VectorBase<Scalar> > > xdot_vec;
  Array<ScalarMag> accuracy_vec;
  interpBuffSource.getPoints(time_vec, &x_vec, &xdot_vec, &accuracy_vec);

  if ( Teuchos::as<int>(this->getVerbLevel()) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
    *out << "Sink buffer range before addPoints = [" << interpBuffSink->getTimeRange().lower() << "," <<
      interpBuffSink->getTimeRange().upper() << "]" << std::endl;
  }

  interpBuffSink->addPoints(time_vec, x_vec, xdot_vec);

  if ( Teuchos::as<int>(this->getVerbLevel()) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
    *out << "Sink buffer range after addPoints = [" << interpBuffSink->getTimeRange().lower() << "," <<
      interpBuffSink->getTimeRange().upper() << "]" << std::endl;
  }

}

//.........这里部分代码省略.........
        // Note that any breakpoint info is also related through this call.
        trialStepCtrlInfo = integrationControlStrategy_->getNextStepControlInfo(
          *stepper_, stepCtrlInfoLast_, currTimeStepIndex_
          );
      }
      else {
        // Take a variable step if we have no control strategy
        trialStepCtrlInfo.stepType = STEP_TYPE_VARIABLE;
        trialStepCtrlInfo.stepSize = NL::max();
      }
    }

    // Print the initial trial step
    if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
      *out << "\nTrial step:\n";
      OSTab tab(out);
      *out << trialStepCtrlInfo;
    }

    // Halt immediately if we where told to do so
    if (trialStepCtrlInfo.stepSize < ST::zero()) {
      if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) )
        *out
          << "\n***"
          << "\n*** NOTICE: The IntegrationControlStrategy object return stepSize < 0.0, halting time integration!"
          << "\n***\n";
      return_val = false;
      break; // Exit the loop immediately!
    }

    // Make sure we don't step past the final time if asked not to
    bool updatedTrialStepCtrlInfo = false;
    {
      const Scalar finalTime = integrationTimeDomain_.upper();
      if (landOnFinalTime_ && trialStepCtrlInfo.stepSize + currStepperTimeRange.upper() > finalTime) {
        if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
          *out << "\nCutting trial step to avoid stepping past final time ...\n";
        trialStepCtrlInfo.stepSize = finalTime - currStepperTimeRange.upper();
        updatedTrialStepCtrlInfo = true;
      }
    }
    
    // Print the modified trial step
    if ( updatedTrialStepCtrlInfo
      && includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) )
    {
      *out << "\nUpdated trial step:\n";
      OSTab tab(out);
      *out << trialStepCtrlInfo;
    }

    //
    // C) Take the step
    //

    // Print step type and size
    if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
      if (trialStepCtrlInfo.stepType == STEP_TYPE_VARIABLE)
        *out << "\nTaking a variable time step with max step size = "
             << trialStepCtrlInfo.stepSize << " ....\n";
      else
        *out << "\nTaking a fixed time step of size = "
             << trialStepCtrlInfo.stepSize << " ....\n";
    }

    // Take step

TEUCHOS_UNIT_TEST( Rythmos_GlobalErrorEstimator, SinCos ) {
  typedef Teuchos::ScalarTraits<double> ST;
  // Forward Solve, storing data in linear interpolation buffer
  int storageLimit = 100;
  double finalTime = 0.1;
  double dt = 0.1;
  RCP<IntegratorBuilder<double> > ib = integratorBuilder<double>();
  {
    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","Backward Euler");
    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<SinCosModel> fwdModel = sinCosModel(true); // implicit formulation
  Thyra::ModelEvaluatorBase::InArgs<double> fwdIC = fwdModel->getNominalValues();
  RCP<Thyra::NonlinearSolverBase<double> > fwdNLSolver = timeStepNonlinearSolver<double>();
  RCP<IntegratorBase<double> > fwdIntegrator = ib->create(fwdModel,fwdIC,fwdNLSolver);
  RCP<const VectorBase<double> > x_final;
  {
    Array<double> time_vec;
    time_vec.push_back(finalTime);
    Array<RCP<const Thyra::VectorBase<double> > > x_final_array;
    fwdIntegrator->getFwdPoints(time_vec,&x_final_array,NULL,NULL);
    x_final = x_final_array[0];
  }
  // Verify x_final is correct
  {
    // Defaults from SinCos Model:
    double f = 1.0;
    double L = 1.0;
    double a = 0.0;
    double x_ic_0 = 0.0;
    double x_ic_1 = 1.0;
    double x_0 = dt/(1.0+std::pow(dt*f/L,2))*(x_ic_0/dt+x_ic_1+dt*std::pow(f/L,2)*a);
    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;
  }
//.........这里部分代码省略.........

TEUCHOS_UNIT_TEST( Rythmos_TimeRange, nonMemberConstructor ) {
  TimeRange<double> tr = timeRange(1.25,3.45);
  TEST_EQUALITY_CONST( tr.isValid(), true );
  TEST_EQUALITY_CONST( tr.lower(), 1.25 );
  TEST_EQUALITY_CONST( tr.upper(), 3.45 );
}