Rythmos - Transient Integration for Differential Equations  Version of the Day
 All Classes Functions Variables Typedefs Pages
List of all members
Rythmos::ForwardSensitivityImplicitModelEvaluator< Scalar > Class Template Reference

Forward sensitivity transient ModelEvaluator subclass. More...

#include <Rythmos_ForwardSensitivityImplicitModelEvaluator.hpp>

Inheritance diagram for Rythmos::ForwardSensitivityImplicitModelEvaluator< Scalar >:
Inheritance graph
[legend]

Constructors/Intializers/Accessors

 ForwardSensitivityImplicitModelEvaluator ()
 
void initializeStructure (const RCP< const Thyra::ModelEvaluator< Scalar > > &stateModel, const int p_index)
 Intialize the with the model structure. More...
 
RCP< const
Thyra::ModelEvaluator< Scalar > > 
getStateModel () const
 
RCP< Thyra::ModelEvaluator
< Scalar > > 
getNonconstStateModel () const
 
int get_p_index () const
 
void setStateIntegrator (const RCP< IntegratorBase< Scalar > > &stateIntegrator, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &stateBasePoint)
 Set the state integrator that will be used to get x and x_dot at various time points. More...
 
void initializePointState (Ptr< StepperBase< Scalar > > stateStepper, bool forceUpToDateW)
 Initialize full state for a single point in time. More...
 

Public functions overridden from ForwardSensitivityModelEvaluatorBase.

Initialize full state for a single point in time.

Parameters
stateBasePoint[in] The base point (x_dot_tilde,x_tilde,t_tilde,...) for which the sensitivities will be computed and represented at time t_tilde. Any sensitivities that are needed must be computed at this same time point. The value of alpha and beta here are ignored.
W_tilde[in,persisting] The composite state derivative matrix computed at the base point stateBasePoint with alpha=coeff_x_dot and beta=coeff_x. This argument can be null, in which case it can be computed internally if needed or not at all.
coeff_x_dot[in] The value of alpha for which W_tilde was computed.
coeff_x[in] The value of beta for which W_tilde was computed.
DfDx_dot[in] The matrix d(f)/d(x_dot) computed at stateBasePoint if available. If this argument is null, then it will be computed internally if needed. The default value is Teuchos::null.
DfDp[in] The matrix d(f)/d(p(p_index)) computed at stateBasePoint if available. If this argument is null, then it will be computed internally if needed. The default value is Teuchos::null.

Preconditions:

  • !is_null(W_tilde)

This function must be called after intializeStructure() and before evalModel() is called. After this function is called, *this model is fully initialized and ready to compute the requested outputs.

void initializeStructureInitCondOnly (const RCP< const Thyra::ModelEvaluator< Scalar > > &stateModel, const RCP< const Thyra::VectorSpaceBase< Scalar > > &p_space)
 
RCP< const
Thyra::DefaultMultiVectorProductVectorSpace
< Scalar > > 
get_s_bar_space () const
 
RCP< const
Thyra::VectorSpaceBase< Scalar > > 
get_p_sens_space () const
 

Public functions overridden from ModelEvaulator.

RCP< const
Thyra::VectorSpaceBase< Scalar > > 
get_x_space () const
 
RCP< const
Thyra::VectorSpaceBase< Scalar > > 
get_f_space () const
 
Thyra::ModelEvaluatorBase::InArgs
< Scalar > 
getNominalValues () const
 
RCP
< Thyra::LinearOpWithSolveBase
< Scalar > > 
create_W () const
 
Thyra::ModelEvaluatorBase::InArgs
< Scalar > 
createInArgs () const
 

Deprecated.

void initializeState (const Thyra::ModelEvaluatorBase::InArgs< Scalar > &stateBasePoint, const RCP< const Thyra::LinearOpWithSolveBase< Scalar > > &W_tilde, const Scalar &coeff_x_dot, const Scalar &coeff_x, const RCP< const Thyra::LinearOpBase< Scalar > > &DfDx_dot=Teuchos::null, const RCP< const Thyra::MultiVectorBase< Scalar > > &DfDp=Teuchos::null)
 Deprecated. More...
 

Additional Inherited Members

Detailed Description

template<class Scalar>
class Rythmos::ForwardSensitivityImplicitModelEvaluator< Scalar >

Forward sensitivity transient ModelEvaluator subclass.

This class provides a very general implemenation of a linear forward sensitivity model evaluator for a DAE.

Introduction

The form of the parameterized state equation is:

  f(x_dot(t),x(t),{p_l},t) = 0, over t = [t0,tf]

  x(t0) = x_init(p)

As shown above, the parameters are assumed to be steady state and can enter through the intial condition and/or through the DAE equation itself.

The forward sensitivity equations written in multi-vector form are:

  F_sens(S_dot,S,t) = d(f)/d(x_dot)*S_dot + d(f)/d(x)*S + d(f)/d(p) = 0, over t = [t0,tf]

  S(t0) = d(x_init)/d(p)

where S is a multi-vector with np columns where each column S(:,j) = d(x)/d(p_j) is the sensitivity of x(t) with respect to the p_j parameter. The sensitivity parameter subvector p here is really just one of the parameter subvectors in the state equation. This index of the parameter subvector for which the sensitivity equations are written for is given by p_index. Note that above d(f)/d(x_dot), d(f)/d(x) and d(f)/d(p are all evaluated at the solution to the state equation (x_dot(t),x(t),t) and are not functions of S_dot or S.

Since the model evaluator interface must be expressed in vector form, the multi-vector form of the forward sensitivity equations is flattened out into:

  f_sens(s_bar_dot(t),s_bar(t),{p_l},t) = 0, over t = [t0,tf]

  s_bar(t0) = s_bar_init

where

  s_bar = [ S(:,0); S(:,0); ...; S(:,np-1) ]

           [ d(f)/d(x_dot)*S_dot(:,0) + d(f)/d(x)*S(:,0) + d(f)/d(p(0))          ]
           [ d(f)/d(x_dot)*S_dot(:,1) + d(f)/d(x)*S(:,1) + d(f)/d(p(1))          ]
  f_sens = [ ...                                                                 ]
           [ d(f)/d(x_dot)*S_dot(:,np-1) + d(f)/d(x)*S(:,np-1) + d(f)/d(p(np-1)) ]

  s_bar_init = [ d(x_init)/d(p(0)); d(x_init)/d(p(1)); ...; d(x_init)/d(p(np-1)) ]

The product vector s_bar is represented as a specialized Thyra::ProductVectorBase subclass object with np "blocks" in terms of a single Thyra::MultiVectorBase object (which has np columns).

Implementation Details

In order to achieve high performance, the representation of the forward sensitivity equations and the computation of the timestep for the sensitivities must reuse as much from the state solve as possible. Here we are most concerned about implicit time stepping methods which compute and use the composite matrix

  W = alpha*d(f)/d(x_dot) + beta*d(f)/d(x)

which is the Jacobian for the nonlinear timestep equation in many methods.

First, let us consider how to represent the forward sensitivity equations using a precomputed LOWSB object

  W_tilde = coeff_x_dot * d(f)/d(x_dot) + coeff_x * d(f)/d(x_dot)

computed at some point (x_dot_tilde,x_tilde,t_tilde,...).

Here is how to evaluate the forward sensitivity equations F_sens using W_tilde:

   F_sens = d(f)/d(x_dot)*S_dot
            + ( (1/coeff_x) * ( W_tilde - coeff_x_dot * d(f)/d(x_dot) ) ) * S
            + d(f)/d(p)

   ==>

   F_sens = (1/coeff_x) * W_tilde * S
            + d(f)/d(x_dot) * ( S_dot - (coeff_x_dot/coeff_x)*S )
            + d(f)/d(p)

Above, we see that in order to compute the multi-vector form of the residual for the forward sensitivity equations F_sens, we must be able to compute something like:

   d(f)/d(x_dot) * V + d(f)/d(p)

This can be done by computing d(f)/d(x_dot) as W with alpha=1.0 and beta=0.0 by calling the underlying state model. Or, a special sensitivity computation could be added to the model evaluator that would compute a generalization of:

   F_sens_var = x_dot_scalar * d(f)/d(x_dot) * V_x_dot
                  + x_scalar * d(f)/d(x) * V_x
                  + p_scalar * d(f)/d(p) * V_p

We could then compute what we need using x_dot_scalar=1.0, x_scalar=0.0, p_scalar=1.0, V_x_dot=V, and V_p=I. For explicit time-stepping methods, this operation could compute the needed sensitivity in one shot. Such an addition to the Thyra::ModelEvaluator interface would be handled through additions to the InArgs and OutArgs classes and the details of which are yet to be worked out.

Up to this point, the only assumption that we have made about the time stepping algorithm is that the sensitivities are only needed and only represented at a single point (x_dot_tilde,x_tilde,t_tilde,...). It is up to the client to ensure that this is indeed the case and it will be the case for many different types of timestepping methods.

In order to get maximum reuse out of a precomputed W_tilde matrix we also have to put in special logic for the evaluation of W_hat with respect to the forward sensitivity residual. When using a precomputed W_tilde from the state timestep computation, the matrix W_hat for the sensitivity residual becomes:

   W_hat = alpha * d(f)/d(x_dot)
           + beta * (1/coeff_x) * ( W_tilde - coeff_x_dot * d(f)/d(x_dot) )

The special logic is that when alpha=coeff_x_dot and beta=coeff_x, then W_hat above simplifies to W_hat=W_tilde as:

   W_hat = coeff_x_dot * d(f)/d(x_dot)
           + coeff_x * (1/coeff_x) * ( W_tilde - coeff_x_dot * d(f)/d(x_dot) )

   ==>

   W_hat = coeff_x_dot * d(f)/d(x_dot) - coeff_x_dot * d(f)/d(x_dot)
           + W_tilde

   ==>

   W_hat = W_tilde

For all of the timestepping methods currently implemented in Rythmos at the time of this writing, alpha=coeff_x_dot and beta=coeff_x will always be true and this is checked for in the code implementation. Because of this, only a single W LOWSB object will be given to a client and only it will be returned. Any other values of alpha and beta requested by the client will thown exceptions. In the future, other values of alpha and beta might be allowed but for now this would only be an error.

Note that there are also a few simplifications in cases where single residual timestepping methods are used. In the case of BDF methods (of which backward Euler is one type), we have:

   S_dot = coeff_x_dot * S_tilde + B_x_dot

   S_dot = coeff_x * S_tilde

In this type of method we see that :

   S_dot - (coeff_x_dot/coeff_x) * S

   ==>

   coeff_x_dot * S_tilde + B_x_dot - (coeff_x_dot/coeff_x) * ( coeff_x * S_tilde )

   ==>

   coeff_x_dot * S_tilde - coeff_x_dot *  S_tilde + B_x_dot

   ==>

   B_x_dot

Therefore, in this type of method, the term involving d(f)/d(x_dot) becomes:

   d(f)/d(x_dot) * ( S_dot - (coeff_x_dot/coeff_x)*S )

   ==>

   d(f)/d(x_dot) * B_x_dot

and is independent of the unknown quantity S_tilde. What this means is that if the residual for the sensitivity equaitions is to be computed multiple times for different values of S_tilde, the term d(f)/d(x_dot) * B_x_dot need only be computed once and can then be reused each time.

ToDo: Finish documention!

Definition at line 302 of file Rythmos_ForwardSensitivityImplicitModelEvaluator.hpp.

Constructor & Destructor Documentation

Member Function Documentation

template<class Scalar >
void Rythmos::ForwardSensitivityImplicitModelEvaluator< Scalar >::initializeStructure ( const RCP< const Thyra::ModelEvaluator< Scalar > > &  stateModel,
const int  p_index 
)
virtual

Intialize the with the model structure.

Parameters
stateModel[in,persisting] The ModelEvaluator that defines the parameterized state model f(x_dot,x,p).
p_index[in] The index of the parameter subvector in stateModel for which sensitivities will be computed for.

This function only intializes the spaces etc. needed to define structure of the problem. *this model object is not fully initialized at this point in that evalModel() will not work yet and will thrown exceptions if called. The function initalizeState() must be called later in order to fully initalize the model.

Implements Rythmos::ForwardSensitivityModelEvaluatorBase< Scalar >.

Definition at line 776 of file Rythmos_ForwardSensitivityImplicitModelEvaluator.hpp.

template<class Scalar >
RCP< const Thyra::ModelEvaluator< Scalar > > Rythmos::ForwardSensitivityImplicitModelEvaluator< Scalar >::getStateModel ( ) const
virtual
template<class Scalar >
RCP< Thyra::ModelEvaluator< Scalar > > Rythmos::ForwardSensitivityImplicitModelEvaluator< Scalar >::getNonconstStateModel ( ) const
virtual
template<class Scalar >
int Rythmos::ForwardSensitivityImplicitModelEvaluator< Scalar >::get_p_index ( ) const
virtual
template<class Scalar >
void Rythmos::ForwardSensitivityImplicitModelEvaluator< Scalar >::setStateIntegrator ( const RCP< IntegratorBase< Scalar > > &  stateIntegrator,
const Thyra::ModelEvaluatorBase::InArgs< Scalar > &  stateBasePoint 
)

Set the state integrator that will be used to get x and x_dot at various time points.

Parameters
stateIntegrator[in,persisting] The integrator that will deliver x and xdot at any valid time t. This integrator must be completely step up and ready to return points from stateIntegrator->getFwdPoints(...).
stateBasePoint[in,persisting] Contains all of the base point information for an evaluation except for t, x, and x_dot. For exmaple, this will contain any parameters or extra input that defines the point.

With a state integrator is set here, the client need not call initializePointState() for individual time points!

Definition at line 620 of file Rythmos_ForwardSensitivityImplicitModelEvaluator.hpp.

template<class Scalar >
void Rythmos::ForwardSensitivityImplicitModelEvaluator< Scalar >::initializePointState ( Ptr< StepperBase< Scalar > >  stateStepper,
bool  forceUpToDateW 
)
virtual

Initialize full state for a single point in time.

Parameters
stateBasePoint[in] The base point (x_dot_tilde,x_tilde,t_tilde,...) for which the sensitivities will be computed and represented at time t_tilde. Any sensitivities that are needed must be computed at this same time point. The value of alpha and beta here are ignored.
W_tilde[in,persisting] The composite state derivative matrix computed at the base point stateBasePoint with alpha=coeff_x_dot and beta=coeff_x. This argument can be null, in which case it can be computed internally if needed or not at all.
coeff_x_dot[in] The value of alpha for which W_tilde was computed.
coeff_x[in] The value of beta for which W_tilde was computed.
DfDx_dot[in] The matrix d(f)/d(x_dot) computed at stateBasePoint if available. If this argument is null, then it will be computed internally if needed. The default value is Teuchos::null.
DfDp[in] The matrix d(f)/d(p(p_index)) computed at stateBasePoint if available. If this argument is null, then it will be computed internally if needed. The default value is Teuchos::null.

Preconditions:

  • nonnull(W_tilde)

This function must be called after intializeStructure() and before evalModel() is called. After this function is called, *this model is fully initialized and ready to compute the requested outputs.

Implements Rythmos::ForwardSensitivityModelEvaluatorBase< Scalar >.

Definition at line 630 of file Rythmos_ForwardSensitivityImplicitModelEvaluator.hpp.

template<class Scalar >
void Rythmos::ForwardSensitivityImplicitModelEvaluator< Scalar >::initializeStructureInitCondOnly ( const RCP< const Thyra::ModelEvaluator< Scalar > > &  stateModel,
const RCP< const Thyra::VectorSpaceBase< Scalar > > &  p_space 
)
virtual
template<class Scalar >
RCP< const Thyra::DefaultMultiVectorProductVectorSpace< Scalar > > Rythmos::ForwardSensitivityImplicitModelEvaluator< Scalar >::get_s_bar_space ( ) const
virtual
template<class Scalar >
RCP< const Thyra::VectorSpaceBase< Scalar > > Rythmos::ForwardSensitivityImplicitModelEvaluator< Scalar >::get_p_sens_space ( ) const
virtual
template<class Scalar >
RCP< const Thyra::VectorSpaceBase< Scalar > > Rythmos::ForwardSensitivityImplicitModelEvaluator< Scalar >::get_x_space ( ) const
template<class Scalar >
RCP< const Thyra::VectorSpaceBase< Scalar > > Rythmos::ForwardSensitivityImplicitModelEvaluator< Scalar >::get_f_space ( ) const
template<class Scalar >
Thyra::ModelEvaluatorBase::InArgs< Scalar > Rythmos::ForwardSensitivityImplicitModelEvaluator< Scalar >::getNominalValues ( ) const
template<class Scalar >
RCP< Thyra::LinearOpWithSolveBase< Scalar > > Rythmos::ForwardSensitivityImplicitModelEvaluator< Scalar >::create_W ( ) const
template<class Scalar >
Thyra::ModelEvaluatorBase::InArgs< Scalar > Rythmos::ForwardSensitivityImplicitModelEvaluator< Scalar >::createInArgs ( ) const
template<class Scalar >
void Rythmos::ForwardSensitivityImplicitModelEvaluator< Scalar >::initializeState ( const Thyra::ModelEvaluatorBase::InArgs< Scalar > &  stateBasePoint,
const RCP< const Thyra::LinearOpWithSolveBase< Scalar > > &  W_tilde,
const Scalar &  coeff_x_dot,
const Scalar &  coeff_x,
const RCP< const Thyra::LinearOpBase< Scalar > > &  DfDx_dot = Teuchos::null,
const RCP< const Thyra::MultiVectorBase< Scalar > > &  DfDp = Teuchos::null 
)

Deprecated.

Definition at line 716 of file Rythmos_ForwardSensitivityImplicitModelEvaluator.hpp.


The documentation for this class was generated from the following file: