9 #ifndef TEMPUS_TEST_VANDERPOL_IMEX_IMPLICITMODEL_IMPL_HPP
10 #define TEMPUS_TEST_VANDERPOL_IMEX_IMPLICITMODEL_IMPL_HPP
12 #include "Teuchos_StandardParameterEntryValidators.hpp"
14 #include "Thyra_DefaultSpmdVectorSpace.hpp"
15 #include "Thyra_DetachedVectorView.hpp"
16 #include "Thyra_DetachedMultiVectorView.hpp"
17 #include "Thyra_DefaultSerialDenseLinearOpWithSolveFactory.hpp"
18 #include "Thyra_DefaultMultiVectorLinearOpWithSolve.hpp"
19 #include "Thyra_DefaultLinearOpSource.hpp"
20 #include "Thyra_VectorStdOps.hpp"
21 #include "Thyra_MultiVectorStdOps.hpp"
22 #include "Thyra_DefaultMultiVectorProductVector.hpp"
26 namespace Tempus_Test {
28 template <
class Scalar>
32 isInitialized_ =
false;
38 acceptModelParams_ =
false;
39 useDfDpAsTangent_ =
false;
47 x_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
48 f_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
50 p_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(np_);
51 dxdp_space_ = Thyra::multiVectorProductVectorSpace(x_space_, np_);
55 setParameterList(pList_);
58 template <
class Scalar>
65 template <
class Scalar>
72 template <
class Scalar>
77 "Error, setupInOutArgs_ must be called first!\n");
78 return nominalValues_;
81 template <
class Scalar>
87 this->get_W_factory();
105 V_V(multivec->col(0).
ptr(), *vec);
111 V_V(multivec->col(1).
ptr(), *vec);
115 Thyra::linearOpWithSolve<Scalar>(*W_factory, matrix);
120 template <
class Scalar>
125 Thyra::createMembers(x_space_, dim_);
129 template <
class Scalar>
134 Thyra::defaultSerialDenseLinearOpWithSolveFactory<Scalar>();
138 template <
class Scalar>
148 template <
class Scalar>
156 template <
class Scalar>
163 using Teuchos::rcp_dynamic_cast;
165 "Error, setupInOutArgs_ must be called first!\n");
168 inArgs.
get_x().assert_not_null();
173 Scalar
eps = epsilon_;
175 if (acceptModelParams_) {
177 if (p_in != Teuchos::null) {
181 if (inArgs.
get_p(1) != Teuchos::null) {
182 dxdp_in = rcp_dynamic_cast<
const DMVPV>(inArgs.
get_p(1),
true)
185 if (inArgs.
get_p(2) != Teuchos::null) {
186 dxdotdp_in = rcp_dynamic_cast<
const DMVPV>(inArgs.
get_p(2),
true)
194 if (acceptModelParams_) {
204 f_out_view[1] = (1.0 - x_in_view[0] * x_in_view[0]) * x_in_view[1] / eps;
208 DfDp_out_view(0, 0) = 0.0;
209 DfDp_out_view(1, 0) =
210 -(1.0 - x_in_view[0] * x_in_view[0]) * x_in_view[1] / (eps * eps);
213 if (useDfDpAsTangent_ && !
is_null(dxdp_in)) {
215 DfDp_out_view(1, 0) +=
216 -2.0 * x_in_view[0] * x_in_view[1] / eps * dxdp(0, 0) +
217 (1.0 - x_in_view[0] * x_in_view[0]) / eps * dxdp(1, 0);
228 -2.0 * beta * x_in_view[0] * x_in_view[1] /
eps;
230 beta * (1.0 - x_in_view[0] * x_in_view[0]) / eps;
237 x_dot_in = inArgs.
get_x_dot().assert_not_null();
243 f_out_view[0] = x_dot_in_view[0];
244 f_out_view[1] = x_dot_in_view[1] -
245 (1.0 - x_in_view[0] * x_in_view[0]) * x_in_view[1] / eps;
249 DfDp_out_view(0, 0) = 0.0;
250 DfDp_out_view(1, 0) =
251 (1.0 - x_in_view[0] * x_in_view[0]) * x_in_view[1] / (eps * eps);
254 if (useDfDpAsTangent_ && !
is_null(dxdotdp_in)) {
256 DfDp_out_view(0, 0) += dxdotdp(0, 0);
257 DfDp_out_view(1, 0) += dxdotdp(1, 0);
259 if (useDfDpAsTangent_ && !
is_null(dxdp_in)) {
261 DfDp_out_view(1, 0) +=
262 2.0 * x_in_view[0] * x_in_view[1] / eps * dxdp(0, 0) -
263 (1.0 - x_in_view[0] * x_in_view[0]) / eps * dxdp(1, 0);
271 W_view(0, 0) = alpha;
274 2.0 * beta * x_in_view[0] * x_in_view[1] /
eps;
275 W_view(1, 1) = alpha - beta * (1.0 - x_in_view[0] * x_in_view[0]) /
282 template <
class Scalar>
286 if (!acceptModelParams_) {
287 return Teuchos::null;
292 else if (l == 1 || l == 2)
294 return Teuchos::null;
297 template <
class Scalar>
301 if (!acceptModelParams_) {
302 return Teuchos::null;
307 p_strings->push_back(
"Model Coefficient: epsilon");
309 p_strings->push_back(
"Model Coefficient: epsilon");
311 p_strings->push_back(
"DxDp");
313 p_strings->push_back(
"Dx_dotDp");
317 template <
class Scalar>
325 template <
class Scalar>
328 if (isInitialized_) {
336 inArgs.
setSupports(Thyra::ModelEvaluatorBase::IN_ARG_x_dot);
337 inArgs.
setSupports(Thyra::ModelEvaluatorBase::IN_ARG_x);
338 inArgs.
setSupports(Thyra::ModelEvaluatorBase::IN_ARG_t);
339 inArgs.
setSupports(Thyra::ModelEvaluatorBase::IN_ARG_alpha);
340 inArgs.
setSupports(Thyra::ModelEvaluatorBase::IN_ARG_beta);
341 if (acceptModelParams_) {
351 outArgs.
setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_f);
352 outArgs.
setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_W_op);
353 if (acceptModelParams_) {
355 outArgs.
setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_DfDp, 0,
356 Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM);
362 nominalValues_ = inArgs_;
364 nominalValues_.set_t(t0_ic_);
366 createMember(x_space_);
369 x_ic_view[0] = x0_ic_;
370 x_ic_view[1] = x1_ic_;
372 nominalValues_.set_x(x_ic);
373 if (acceptModelParams_) {
375 createMember(p_space_);
378 p_ic_view[0] = epsilon_;
380 nominalValues_.set_p(0, p_ic);
383 createMember(x_space_);
386 x_dot_ic_view[0] = x1_ic_;
387 x_dot_ic_view[1] = ((1.0 - x0_ic_ * x0_ic_) * x1_ic_ - x0_ic_) / epsilon_;
389 nominalValues_.set_x_dot(x_dot_ic);
392 isInitialized_ =
true;
395 template <
class Scalar>
403 if (paramList != Teuchos::null) tmpPL = paramList;
405 this->setMyParamList(tmpPL);
407 bool acceptModelParams = get<bool>(*pl,
"Accept model parameters");
408 bool haveIC = get<bool>(*pl,
"Provide nominal values");
409 bool useDfDpAsTangent = get<bool>(*pl,
"Use DfDp as Tangent");
410 if ((acceptModelParams != acceptModelParams_) || (haveIC != haveIC_)) {
411 isInitialized_ =
false;
413 acceptModelParams_ = acceptModelParams;
415 useDfDpAsTangent_ = useDfDpAsTangent;
416 epsilon_ = get<Scalar>(*pl,
"Coeff epsilon");
417 x0_ic_ = get<Scalar>(*pl,
"IC x0");
418 x1_ic_ = get<Scalar>(*pl,
"IC x1");
419 t0_ic_ = get<Scalar>(*pl,
"IC t0");
423 template <
class Scalar>
430 pl->
set(
"Accept model parameters",
false);
431 pl->
set(
"Provide nominal values",
true);
432 pl->
set(
"Use DfDp as Tangent",
false);
433 Teuchos::setDoubleParameter(
"Coeff epsilon", 1.0e-06,
434 "Coefficient a in model", &*pl);
435 Teuchos::setDoubleParameter(
"IC x0", 2.0,
"Initial Condition for x0", &*pl);
436 Teuchos::setDoubleParameter(
"IC x1", 0.0,
"Initial Condition for x1", &*pl);
437 Teuchos::setDoubleParameter(
"IC t0", 0.0,
"Initial time t0", &*pl);
444 #endif // TEMPUS_TEST_VANDERPOL_IMEX_IMPLICITMODEL_IMPL_HPP
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int l) const
void setupInOutArgs_() const
Teuchos::RCP< Thyra::LinearOpWithSolveBase< Scalar > > create_W() const
bool is_null(const boost::shared_ptr< T > &p)
Derivative< Scalar > get_DfDp(int l) const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
RCP< const VectorBase< Scalar > > get_x_dot() const
void setModelEvalDescription(const std::string &modelEvalDescription)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
VanDerPol_IMEX_ImplicitModel(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
void set_Np_Ng(int Np, int Ng)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
Evaluation< VectorBase< Scalar > > get_f() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_g_space(int j) const
#define TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE(index, lower_inclusive, upper_exclusive)
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
void setSupports(EInArgsMembers arg, bool supports=true)
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs_bar, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs_bar) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
void setSupports(EOutArgsMembers arg, bool supports=true)
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const ¶mList)
void setModelEvalDescription(const std::string &modelEvalDescription)
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
RCP< LinearOpBase< Scalar > > get_W_op() const
RCP< const VectorBase< Scalar > > get_x() const
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
RCP< const VectorBase< Scalar > > get_p(int l) const
RCP< MultiVectorBase< Scalar > > getMultiVector() const