Skip to content

Commit

Permalink
Merge Pull Request #8814 from trilinos/Trilinos/dirk-dahlquist
Browse files Browse the repository at this point in the history
Automatically Merged using Trilinos Pull Request AutoTester
PR Title: Dirk dahlquist
PR Author: sconde
  • Loading branch information
trilinos-autotester authored Mar 15, 2021
2 parents 09f7d11 + 29f6d76 commit 5434315
Show file tree
Hide file tree
Showing 5 changed files with 502 additions and 25 deletions.
4 changes: 2 additions & 2 deletions packages/tempus/src/Tempus_StepperRKButcherTableau.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -4748,8 +4748,8 @@ createStepperSDIRK_5Stage5thOrder(
* & b^{*T}
* \end{array}
* \;\;\;\;\mbox{ where }\;\;\;\;
* \begin{array}{c|cccc} 0 & 0 & \\
* 1 & -1 & 1 \\ \hline
* \begin{array}{c|cccc} 1 & 1 & \\
* 0 & -1 & 1 \\ \hline
* & 1/2 & 1/2 \\
* & 1 & 0 \end{array}
* \f]
Expand Down
25 changes: 16 additions & 9 deletions packages/tempus/test/TestModels/DahlquistTestModel_decl.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -58,14 +58,14 @@ class DahlquistTestModel
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;
//Teuchos::RCP<Thyra::LinearOpBase<Scalar> > create_W_op() const;
//Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > get_W_factory() const;
Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> > create_W() const;
Teuchos::RCP<Thyra::LinearOpBase<Scalar> > create_W_op() const;
Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > get_W_factory() const;
Thyra::ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;

//Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > get_p_space(int l) const;
//Teuchos::RCP<const Teuchos::Array<std::string> > get_p_names(int l) const;
//Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > get_g_space(int j) const;
Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > get_p_space(int l) const;
Teuchos::RCP<const Teuchos::Array<std::string> > get_p_names(int l) const;
Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > get_g_space(int j) const;
//@}

private:
Expand All @@ -80,6 +80,13 @@ class DahlquistTestModel
//@}

private:
int dim_;
int Np_;
int np_;
int Ng_;
int ng_;
bool haveIC_;
bool acceptModelParams_;
Scalar lambda_;
bool includeXDot_;

Expand All @@ -89,9 +96,9 @@ class DahlquistTestModel
mutable Thyra::ModelEvaluatorBase::InArgs<Scalar> nominalValues_;
Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > x_space_;
Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > f_space_;
//Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > p_space_;
//Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > g_space_;
//Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > DxDp_space_;
Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > p_space_;
Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > g_space_;
Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > DxDp_space_;

Scalar xIC_; ///< Initial condition for x.
Scalar xDotIC_; ///< Initial condition for xDot.
Expand Down
175 changes: 163 additions & 12 deletions packages/tempus/test/TestModels/DahlquistTestModel_impl.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -11,13 +11,13 @@

#include "Thyra_DefaultSpmdVectorSpace.hpp"
#include "Thyra_DetachedVectorView.hpp"
//#include "Thyra_DetachedMultiVectorView.hpp"
//#include "Thyra_DefaultSerialDenseLinearOpWithSolveFactory.hpp"
//#include "Thyra_DefaultMultiVectorLinearOpWithSolve.hpp"
//#include "Thyra_DefaultLinearOpSource.hpp"
//#include "Thyra_VectorStdOps.hpp"
//#include "Thyra_MultiVectorStdOps.hpp"
//#include "Thyra_DefaultMultiVectorProductVector.hpp"
#include "Thyra_DetachedMultiVectorView.hpp"
#include "Thyra_DefaultSerialDenseLinearOpWithSolveFactory.hpp"
#include "Thyra_DefaultMultiVectorLinearOpWithSolve.hpp"
#include "Thyra_DefaultLinearOpSource.hpp"
#include "Thyra_VectorStdOps.hpp"
#include "Thyra_MultiVectorStdOps.hpp"
#include "Thyra_DefaultMultiVectorProductVector.hpp"

//#include <iostream>

Expand Down Expand Up @@ -45,11 +45,17 @@ constructDahlquistTestModel(Scalar lambda, bool includeXDot)
isInitialized_ = false;
xIC_ = Scalar( 1.0);
xDotIC_ = Scalar(lambda);
int dim = 1;
dim_ = 1;
haveIC_ = true;
Np_ = 1;
np_ = 1;
Ng_ = 1;
ng_ = dim_;
acceptModelParams_ = false;

// Create x_space and f_space
x_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim);
f_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim);
x_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
f_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);

using Teuchos::RCP;
typedef Thyra::ModelEvaluatorBase MEB;
Expand All @@ -60,6 +66,11 @@ constructDahlquistTestModel(Scalar lambda, bool includeXDot)
inArgs.setSupports( MEB::IN_ARG_t );
inArgs.setSupports( MEB::IN_ARG_x );
inArgs.setSupports( MEB::IN_ARG_x_dot );

inArgs.setSupports( MEB::IN_ARG_beta );
inArgs.setSupports( MEB::IN_ARG_alpha );
if (acceptModelParams_) inArgs.set_Np( Np_ );

inArgs_ = inArgs;
}

Expand All @@ -68,6 +79,17 @@ constructDahlquistTestModel(Scalar lambda, bool includeXDot)
MEB::OutArgsSetup<Scalar> outArgs;
outArgs.setModelEvalDescription(this->description());
outArgs.setSupports( MEB::OUT_ARG_f );

outArgs.setSupports( MEB::OUT_ARG_W_op );
if (acceptModelParams_) {
outArgs.set_Np_Ng(Np_,Ng_);
outArgs.setSupports( MEB::OUT_ARG_DfDp,0,
MEB::DERIV_MV_JACOBIAN_FORM );
outArgs.setSupports( MEB::OUT_ARG_DgDp,0,0,
MEB::DERIV_MV_JACOBIAN_FORM );
outArgs.setSupports( MEB::OUT_ARG_DgDx,0,
MEB::DERIV_MV_GRADIENT_FORM );
}
outArgs_ = outArgs;
}

Expand Down Expand Up @@ -180,12 +202,18 @@ evalModelImpl(
const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
) const
{
typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
using Teuchos::RCP;
using Thyra::VectorBase;
using Thyra::MultiVectorBase;
using Teuchos::rcp_dynamic_cast;
TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
"Error, setupInOutArgs_ must be called first!\n");

const RCP<const VectorBase<Scalar> > x_in = inArgs.get_x().assert_not_null();
Thyra::ConstDetachedVectorView<Scalar> x_in_view( *x_in );
const RCP<VectorBase<Scalar> > f_out = outArgs.get_f();
const RCP<Thyra::LinearOpBase<Scalar> > W_out = outArgs.get_W_op();

if (inArgs.get_x_dot().is_null()) {

Expand All @@ -197,13 +225,136 @@ evalModelImpl(
TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error,
"Error -- Dahlquist Test Model requires f_out!\n");
}

if (!is_null(W_out)) {
RCP<Thyra::MultiVectorBase<Scalar> > matrix =
Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,true);
Thyra::DetachedMultiVectorView<Scalar> matrix_view( *matrix );
matrix_view(0,0) = lambda_;
}

} else {
TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error,
"Error -- Dahlquist Test Model only setup for explicit ODE!\n");

// Evaluate the implicit ODE f(xdot, x, t) [=0]
RCP<const VectorBase<Scalar> > x_dot_in;
x_dot_in = inArgs.get_x_dot().assert_not_null();
Scalar alpha = inArgs.get_alpha();
Scalar beta = inArgs.get_beta();

if (!is_null(f_out)) {
Thyra::DetachedVectorView<Scalar> f_out_view( *f_out );
Thyra::ConstDetachedVectorView<Scalar> x_dot_in_view( *x_dot_in );
f_out_view[0] = x_dot_in_view[0] - lambda_*x_in_view[0];
}

if (!is_null(W_out)) {
RCP<Thyra::MultiVectorBase<Scalar> > matrix =
Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(W_out,true);
Thyra::DetachedMultiVectorView<Scalar> matrix_view( *matrix );
matrix_view(0,0) = alpha - beta*lambda_; // d(f0)/d(x0_n)
// Note: alpha = d(xdot)/d(x_n) and beta = d(x)/d(x_n)
}

}

}

template<class Scalar>
Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
DahlquistTestModel<Scalar>::
create_W() const
{
using Teuchos::RCP;
RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > W_factory = this->get_W_factory();
RCP<Thyra::LinearOpBase<Scalar> > matrix = this->create_W_op();
{
RCP<Thyra::MultiVectorBase<Scalar> > multivec = Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(matrix,true);
{
RCP<Thyra::VectorBase<Scalar> > vec = Thyra::createMember(x_space_);
{
Thyra::DetachedVectorView<Scalar> vec_view( *vec );
vec_view[0] = lambda_;
}
V_V(multivec->col(0).ptr(),*vec);
}
}
RCP<Thyra::LinearOpWithSolveBase<Scalar> > W =
Thyra::linearOpWithSolve<Scalar>(
*W_factory,
matrix
);
return W;
}

template<class Scalar>
Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
DahlquistTestModel<Scalar>::
create_W_op() const
{
//const int dim_ = 1;
Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > matrix = Thyra::createMembers(x_space_, dim_);
return(matrix);
}


template<class Scalar>
Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
DahlquistTestModel<Scalar>::
get_W_factory() const
{
Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > W_factory =
Thyra::defaultSerialDenseLinearOpWithSolveFactory<Scalar>();
return W_factory;
}

template<class Scalar>
Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
DahlquistTestModel<Scalar>::
get_p_space(int l) const
{
if (!acceptModelParams_) {
return Teuchos::null;
}
TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
if (l == 0)
return p_space_;
else if (l == 1 || l == 2)
return DxDp_space_;
return Teuchos::null;
}


template<class Scalar>
Teuchos::RCP<const Teuchos::Array<std::string> >
DahlquistTestModel<Scalar>::
get_p_names(int l) const
{
if (!acceptModelParams_) {
return Teuchos::null;
}
TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
Teuchos::RCP<Teuchos::Array<std::string> > p_strings =
Teuchos::rcp(new Teuchos::Array<std::string>());
if (l == 0) {
p_strings->push_back("Model Coefficient: a");
//p_strings->push_back("Model Coefficient: f");
//p_strings->push_back("Model Coefficient: L");
}
else if (l == 1)
p_strings->push_back("DxDp");
else if (l == 2)
p_strings->push_back("Dx_dotDp");
return p_strings;
}

template<class Scalar>
Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
DahlquistTestModel<Scalar>::
get_g_space(int j) const
{
TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( j, 0, Ng_ );
return g_space_;
}

} // namespace Tempus_Test
#endif // TEMPUS_TEST_DAHLQUIST_TEST_MODEL_IMPL_HPP
Loading

0 comments on commit 5434315

Please sign in to comment.