Rythmos::ForwardSensitivityModelEvaluator< Scalar > Class Template Reference

Forward sensitivity transient ModelEvaluator subclass. More...

#include <Rythmos_ForwardSensitivityModelEvaluator.hpp>

Inheritance diagram for Rythmos::ForwardSensitivityModelEvaluator< Scalar >:

[legend]
List of all members.

Constructors/Intializers/Accessors

 ForwardSensitivityModelEvaluator ()
 
void initializeStructure (const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &stateModel, const int p_index)
 Intialize the with the model structure.
Teuchos::RCP< const Thyra::ModelEvaluator<
Scalar > > 
getStateModel () const
 
int get_p_index () const
 
void initializeState (const Thyra::ModelEvaluatorBase::InArgs< Scalar > &stateBasePoint, const Teuchos::RCP< const Thyra::LinearOpWithSolveBase< Scalar > > &W_tilde, const Scalar &coeff_x_dot, const Scalar &coeff_x, const Teuchos::RCP< const Thyra::LinearOpBase< Scalar > > &DfDx_dot=Teuchos::null, const Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > &DfDp=Teuchos::null)
 Initialize full state.

Public functions overridden from ModelEvaulator.

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

Detailed Description

template<class Scalar>
class Rythmos::ForwardSensitivityModelEvaluator< 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 types 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 297 of file Rythmos_ForwardSensitivityModelEvaluator.hpp.


Constructor & Destructor Documentation

template<class Scalar>
Rythmos::ForwardSensitivityModelEvaluator< Scalar >::ForwardSensitivityModelEvaluator (  ) 

Definition at line 465 of file Rythmos_ForwardSensitivityModelEvaluator.hpp.


Member Function Documentation

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

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.

Definition at line 471 of file Rythmos_ForwardSensitivityModelEvaluator.hpp.

template<class Scalar>
Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > Rythmos::ForwardSensitivityModelEvaluator< Scalar >::getStateModel (  )  const

Definition at line 516 of file Rythmos_ForwardSensitivityModelEvaluator.hpp.

template<class Scalar>
int Rythmos::ForwardSensitivityModelEvaluator< Scalar >::get_p_index (  )  const

Definition at line 523 of file Rythmos_ForwardSensitivityModelEvaluator.hpp.

template<class Scalar>
void Rythmos::ForwardSensitivityModelEvaluator< Scalar >::initializeState ( const Thyra::ModelEvaluatorBase::InArgs< Scalar > &  stateBasePoint,
const Teuchos::RCP< const Thyra::LinearOpWithSolveBase< Scalar > > &  W_tilde,
const Scalar &  coeff_x_dot,
const Scalar &  coeff_x,
const Teuchos::RCP< const Thyra::LinearOpBase< Scalar > > &  DfDx_dot = Teuchos::null,
const Teuchos::RCP< const Thyra::MultiVectorBase< Scalar > > &  DfDp = Teuchos::null 
)

Initialize full state.

Parameters:
stateBasePoint [in] The base point (x_dot_tilde,x_tilde,t_tilde,...) for which the sensitivities will be computed and represented at. Any sensitivities that are needed are computed at this 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 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 value is Teuchos::null.
Preconditions:

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.

Definition at line 530 of file Rythmos_ForwardSensitivityModelEvaluator.hpp.

template<class Scalar>
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > Rythmos::ForwardSensitivityModelEvaluator< Scalar >::get_x_space (  )  const [virtual]

Implements Thyra::ModelEvaluator< Scalar >.

Definition at line 588 of file Rythmos_ForwardSensitivityModelEvaluator.hpp.

template<class Scalar>
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > Rythmos::ForwardSensitivityModelEvaluator< Scalar >::get_f_space (  )  const [virtual]

Implements Thyra::ModelEvaluator< Scalar >.

Definition at line 596 of file Rythmos_ForwardSensitivityModelEvaluator.hpp.

template<class Scalar>
Thyra::ModelEvaluatorBase::InArgs< Scalar > Rythmos::ForwardSensitivityModelEvaluator< Scalar >::getNominalValues (  )  const [virtual]

Reimplemented from Thyra::StateFuncModelEvaluatorBase< Scalar >.

Definition at line 604 of file Rythmos_ForwardSensitivityModelEvaluator.hpp.

template<class Scalar>
Teuchos::RCP< Thyra::LinearOpWithSolveBase< Scalar > > Rythmos::ForwardSensitivityModelEvaluator< Scalar >::create_W (  )  const [virtual]

Reimplemented from Thyra::StateFuncModelEvaluatorBase< Scalar >.

Definition at line 612 of file Rythmos_ForwardSensitivityModelEvaluator.hpp.

template<class Scalar>
Thyra::ModelEvaluatorBase::InArgs< Scalar > Rythmos::ForwardSensitivityModelEvaluator< Scalar >::createInArgs (  )  const [virtual]

Implements Thyra::ModelEvaluator< Scalar >.

Definition at line 620 of file Rythmos_ForwardSensitivityModelEvaluator.hpp.


The documentation for this class was generated from the following file:
Generated on Tue Oct 20 12:46:09 2009 for Rythmos - Transient Integration for Differential Equations by doxygen 1.4.7