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

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() );
}

Scalar translateTimeRange(
    Scalar t,
    const TimeRange<Scalar>& sourceRange,
    const TimeRange<Scalar>& destinationRange
    ) {
  Scalar r = destinationRange.length()/sourceRange.length();
  return r*t+destinationRange.lower()-r*sourceRange.lower();
}

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::assertNoTimePointsBeforeCurrentTimeRange(
  const InterpolationBufferBase<Scalar> &interpBuffer,
  const Array<Scalar>& time_vec,
  const int &startingTimePointIndex
  )
{
  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(
        time_vec[i] < currentTimeRange.lower(), std::out_of_range,
        "Error, time_vec["<<i<<"] = " << time_vec[i] << " < currentTimeRange.lower() = "
        << currentTimeRange.lower() << " for " << interpBuffer.description() << "!"
        );
    }
  }
}

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() << "]!"
        );
    }
  }
}

bool Rythmos::getCurrentPoints(
  const InterpolationBufferBase<Scalar> &interpBuffer,
  const Array<Scalar>& time_vec,
  Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
  Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
  int *nextTimePointIndex_inout
  )
{

  typedef ScalarTraits<Scalar> ST;
  using Teuchos::as;

  const int numTotalTimePoints = time_vec.size();

  // Validate input
#ifdef RYTHMOS_DEBUG
  TEST_FOR_EXCEPT(nextTimePointIndex_inout==0);
  TEUCHOS_ASSERT( 0 <= *nextTimePointIndex_inout && *nextTimePointIndex_inout < numTotalTimePoints );
  TEUCHOS_ASSERT( x_vec == 0 || as<int>(x_vec->size()) == numTotalTimePoints );
  TEUCHOS_ASSERT( xdot_vec == 0 || as<int>(xdot_vec->size()) == numTotalTimePoints );
#endif // RYTHMOS_DEBUG

  int &nextTimePointIndex = *nextTimePointIndex_inout;
  const int initNextTimePointIndex = nextTimePointIndex;

  const TimeRange<Scalar> currentTimeRange = interpBuffer.getTimeRange();
  
  if (currentTimeRange.length() >= ST::zero()) {

    // Load a temp array with all of the current time points that fall in the
    // current time range.
    Array<Scalar> current_time_vec;
    { // scope for i to remove shadow warning.
      int i;
      for ( i = 0; i < numTotalTimePoints-nextTimePointIndex; ++i ) {
        const Scalar t = time_vec[nextTimePointIndex];
#ifdef RYTHMOS_DEBUG
        TEUCHOS_ASSERT( t >= currentTimeRange.lower() );
#endif // RYTHMOS_DEBUG
        if ( currentTimeRange.isInRange(t) ) {
          ++nextTimePointIndex;
          current_time_vec.push_back(t);
        }
        else {
          break;
        }
      }
#ifdef RYTHMOS_DEBUG
      // Here I am just checking that the loop worked as expected with the data
      // in the current time range all comming first.
      TEUCHOS_ASSERT( nextTimePointIndex-initNextTimePointIndex == i );
#endif
    }

    // Get points in current time range if any such points exist

    const int numCurrentTimePoints = current_time_vec.size();

    if ( numCurrentTimePoints > 0 ) {

      // Get the state(s) for current time points from the stepper and put
      // them into temp arrays
      Array<RCP<const Thyra::VectorBase<Scalar> > > current_x_vec;
      Array<RCP<const Thyra::VectorBase<Scalar> > > current_xdot_vec;
      if (x_vec || xdot_vec) {
        interpBuffer.getPoints(
          current_time_vec,
          x_vec ? &current_x_vec : 0,
          xdot_vec ? &current_xdot_vec : 0,
          0 // accuracy_vec
          );
      }

      // Copy the gotten x and xdot vectors from the temp arrays to the output
      // arrays.
      for ( int i = initNextTimePointIndex; i < nextTimePointIndex; ++i ) {
        if (x_vec)
          (*x_vec)[i] = current_x_vec[i-initNextTimePointIndex];
        if (xdot_vec)
          (*xdot_vec)[i] = current_xdot_vec[i-initNextTimePointIndex];
      }

    }

  }

  return ( nextTimePointIndex == initNextTimePointIndex ? false : true );

}

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;
  }

}

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 );
}