ModelEvaluator subclass.
More...
#include <Rythmos_ForwardSensitivityModelEvaluator.hpp>
Inheritance diagram for Rythmos::ForwardSensitivityModelEvaluator< Scalar >:
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 |
| | |
ModelEvaluator subclass.
This class provides a very general implemenation of a linear forward sensitivity model evaluator for a DAE.
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).
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.
| Rythmos::ForwardSensitivityModelEvaluator< Scalar >::ForwardSensitivityModelEvaluator | ( | ) |
| void Rythmos::ForwardSensitivityModelEvaluator< Scalar >::initializeStructure | ( | const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > & | stateModel, | |
| const int | p_index | |||
| ) |
Intialize the with the model structure.
| 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 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.
| Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > Rythmos::ForwardSensitivityModelEvaluator< Scalar >::getStateModel | ( | ) | const |
| int Rythmos::ForwardSensitivityModelEvaluator< Scalar >::get_p_index | ( | ) | const |
| 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.
| 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. |
!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.
Definition at line 530 of file Rythmos_ForwardSensitivityModelEvaluator.hpp.
| 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.
| 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.
| Thyra::ModelEvaluatorBase::InArgs< Scalar > Rythmos::ForwardSensitivityModelEvaluator< Scalar >::getNominalValues | ( | ) | const [virtual] |
Reimplemented from Thyra::StateFuncModelEvaluatorBase< Scalar >.
Definition at line 604 of file Rythmos_ForwardSensitivityModelEvaluator.hpp.
| 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.
| Thyra::ModelEvaluatorBase::InArgs< Scalar > Rythmos::ForwardSensitivityModelEvaluator< Scalar >::createInArgs | ( | ) | const [virtual] |
Implements Thyra::ModelEvaluator< Scalar >.
Definition at line 620 of file Rythmos_ForwardSensitivityModelEvaluator.hpp.
1.4.7