From 22e415734d98624112456f97b7641ba3bcdaa349 Mon Sep 17 00:00:00 2001 From: Sidafa Conde Date: Tue, 23 Feb 2021 14:57:08 -0500 Subject: [PATCH 1/2] Tempus: TestModels - update DahlquistTestModel for DIRK stepper - define implicit ME operators - unit_test - DIRK - define the implicit residual computation - group switch cases to reduce code duplication - DIRK add ESDIRK_Trap to the test --- .../src/Tempus_StepperRKButcherTableau.hpp | 4 +- .../TestModels/DahlquistTestModel_decl.hpp | 25 +- .../TestModels/DahlquistTestModel_impl.hpp | 175 ++++++++- packages/tempus/unit_test/CMakeLists.txt | 7 + .../tempus/unit_test/Tempus_UnitTest_DIRK.cpp | 343 ++++++++++++++++++ 5 files changed, 531 insertions(+), 23 deletions(-) create mode 100644 packages/tempus/unit_test/Tempus_UnitTest_DIRK.cpp diff --git a/packages/tempus/src/Tempus_StepperRKButcherTableau.hpp b/packages/tempus/src/Tempus_StepperRKButcherTableau.hpp index efe2747f31da..744822252e7d 100644 --- a/packages/tempus/src/Tempus_StepperRKButcherTableau.hpp +++ b/packages/tempus/src/Tempus_StepperRKButcherTableau.hpp @@ -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] diff --git a/packages/tempus/test/TestModels/DahlquistTestModel_decl.hpp b/packages/tempus/test/TestModels/DahlquistTestModel_decl.hpp index f13ca0806491..dad4e9bbff9f 100644 --- a/packages/tempus/test/TestModels/DahlquistTestModel_decl.hpp +++ b/packages/tempus/test/TestModels/DahlquistTestModel_decl.hpp @@ -58,14 +58,14 @@ class DahlquistTestModel Teuchos::RCP > get_x_space() const; Teuchos::RCP > get_f_space() const; Thyra::ModelEvaluatorBase::InArgs getNominalValues() const; - //Teuchos::RCP > create_W() const; - //Teuchos::RCP > create_W_op() const; - //Teuchos::RCP > get_W_factory() const; + Teuchos::RCP > create_W() const; + Teuchos::RCP > create_W_op() const; + Teuchos::RCP > get_W_factory() const; Thyra::ModelEvaluatorBase::InArgs createInArgs() const; - //Teuchos::RCP > get_p_space(int l) const; - //Teuchos::RCP > get_p_names(int l) const; - //Teuchos::RCP > get_g_space(int j) const; + Teuchos::RCP > get_p_space(int l) const; + Teuchos::RCP > get_p_names(int l) const; + Teuchos::RCP > get_g_space(int j) const; //@} private: @@ -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_; @@ -89,9 +96,9 @@ class DahlquistTestModel mutable Thyra::ModelEvaluatorBase::InArgs nominalValues_; Teuchos::RCP > x_space_; Teuchos::RCP > f_space_; - //Teuchos::RCP > p_space_; - //Teuchos::RCP > g_space_; - //Teuchos::RCP > DxDp_space_; + Teuchos::RCP > p_space_; + Teuchos::RCP > g_space_; + Teuchos::RCP > DxDp_space_; Scalar xIC_; ///< Initial condition for x. Scalar xDotIC_; ///< Initial condition for xDot. diff --git a/packages/tempus/test/TestModels/DahlquistTestModel_impl.hpp b/packages/tempus/test/TestModels/DahlquistTestModel_impl.hpp index 10e010734e55..dff7766aa4e6 100644 --- a/packages/tempus/test/TestModels/DahlquistTestModel_impl.hpp +++ b/packages/tempus/test/TestModels/DahlquistTestModel_impl.hpp @@ -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 @@ -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(dim); - f_space_ = Thyra::defaultSpmdVectorSpace(dim); + x_space_ = Thyra::defaultSpmdVectorSpace(dim_); + f_space_ = Thyra::defaultSpmdVectorSpace(dim_); using Teuchos::RCP; typedef Thyra::ModelEvaluatorBase MEB; @@ -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; } @@ -68,6 +79,17 @@ constructDahlquistTestModel(Scalar lambda, bool includeXDot) MEB::OutArgsSetup 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; } @@ -180,12 +202,18 @@ evalModelImpl( const Thyra::ModelEvaluatorBase::OutArgs &outArgs ) const { + typedef Thyra::DefaultMultiVectorProductVector 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 > x_in = inArgs.get_x().assert_not_null(); Thyra::ConstDetachedVectorView x_in_view( *x_in ); const RCP > f_out = outArgs.get_f(); + const RCP > W_out = outArgs.get_W_op(); if (inArgs.get_x_dot().is_null()) { @@ -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 > matrix = + Teuchos::rcp_dynamic_cast >(W_out,true); + Thyra::DetachedMultiVectorView 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 > 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 f_out_view( *f_out ); + Thyra::ConstDetachedVectorView 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 > matrix = + Teuchos::rcp_dynamic_cast >(W_out,true); + Thyra::DetachedMultiVectorView 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 +Teuchos::RCP > +DahlquistTestModel:: +create_W() const +{ + using Teuchos::RCP; + RCP > W_factory = this->get_W_factory(); + RCP > matrix = this->create_W_op(); + { + RCP > multivec = Teuchos::rcp_dynamic_cast >(matrix,true); + { + RCP > vec = Thyra::createMember(x_space_); + { + Thyra::DetachedVectorView vec_view( *vec ); + vec_view[0] = lambda_; + } + V_V(multivec->col(0).ptr(),*vec); + } + } + RCP > W = + Thyra::linearOpWithSolve( + *W_factory, + matrix + ); + return W; +} + +template +Teuchos::RCP > +DahlquistTestModel:: +create_W_op() const +{ + //const int dim_ = 1; + Teuchos::RCP > matrix = Thyra::createMembers(x_space_, dim_); + return(matrix); +} + + +template +Teuchos::RCP > +DahlquistTestModel:: +get_W_factory() const +{ + Teuchos::RCP > W_factory = + Thyra::defaultSerialDenseLinearOpWithSolveFactory(); + return W_factory; +} + +template +Teuchos::RCP > +DahlquistTestModel:: +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 +Teuchos::RCP > +DahlquistTestModel:: +get_p_names(int l) const +{ + if (!acceptModelParams_) { + return Teuchos::null; + } + TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ ); + Teuchos::RCP > p_strings = + Teuchos::rcp(new Teuchos::Array()); + 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 +Teuchos::RCP > +DahlquistTestModel:: +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 diff --git a/packages/tempus/unit_test/CMakeLists.txt b/packages/tempus/unit_test/CMakeLists.txt index 8d79cca20684..931aa47ed89a 100644 --- a/packages/tempus/unit_test/CMakeLists.txt +++ b/packages/tempus/unit_test/CMakeLists.txt @@ -414,6 +414,13 @@ TRIBITS_ADD_EXECUTABLE_AND_TEST( NUM_MPI_PROCS 1 ) +TRIBITS_ADD_EXECUTABLE_AND_TEST( + UnitTest_DIRK + SOURCES Tempus_UnitTest_DIRK.cpp ${TEUCHOS_STD_UNIT_TEST_MAIN} + TESTONLYLIBS tempus_test_models + NUM_MPI_PROCS 1 + ) + TRIBITS_COPY_FILES_TO_BINARY_DIR(UnitTest_NewmarkImplicitAForm_CopyFiles DEST_FILES Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder.xml Tempus_NewmarkExplicitAForm_HarmonicOscillator_Damped.xml EXEDEPS UnitTest_NewmarkImplicitAForm diff --git a/packages/tempus/unit_test/Tempus_UnitTest_DIRK.cpp b/packages/tempus/unit_test/Tempus_UnitTest_DIRK.cpp new file mode 100644 index 000000000000..352cc55a94e3 --- /dev/null +++ b/packages/tempus/unit_test/Tempus_UnitTest_DIRK.cpp @@ -0,0 +1,343 @@ +// @HEADER +// **************************************************************************** +// Tempus: Copyright (2017) Sandia Corporation +// +// Distributed under BSD 3-clause license (See accompanying file Copyright.txt) +// **************************************************************************** +// @HEADER + +#include "Teuchos_UnitTestHarness.hpp" + +#include "../TestModels/DahlquistTestModel.hpp" + +#include "Tempus_SolutionHistory.hpp" +#include "Tempus_StepperRKButcherTableau.hpp" +#include "Tempus_StepperRKBase.hpp" +#include "Tempus_StepperRKModifierBase.hpp" +#include "Tempus_StepperRKModifierXBase.hpp" + + +namespace Tempus_Unit_Test { + +using Teuchos::RCP; +using Teuchos::rcp; +using Teuchos::rcp_const_cast; +using Teuchos::rcp_dynamic_cast; + + +class StepperRKModifierSDIRK21Test + : virtual public Tempus::StepperRKModifierBase +{ +public: + + /// Constructor + StepperRKModifierSDIRK21Test(Teuchos::FancyOStream &Out, bool &Success) + : out(Out), success(Success) + {} + + /// Destructor + virtual ~StepperRKModifierSDIRK21Test(){} + + /// Test the modify RK Stepper at the Action Locations. + virtual void modify( + Teuchos::RCP > sh, + Teuchos::RCP > stepper, + const typename Tempus::StepperRKAppAction::ACTION_LOCATION actLoc) + { + const double relTol = 1.0e-14; + auto stageNumber = stepper->getStageNumber(); + Teuchos::SerialDenseVector c = stepper->getTableau()->c(); + + auto currentState = sh->getCurrentState(); + auto workingState = sh->getWorkingState(); + const double dt = workingState->getTimeStep(); + double time = currentState->getTime(); + if (stageNumber >= 0) time += c(stageNumber)*dt; + + auto x = workingState->getX(); + auto xDot = workingState->getXDot(); + if (xDot == Teuchos::null) xDot = stepper->getStepperXDot(); + + switch (actLoc) { + case StepperRKAppAction::BEGIN_STEP: + { + { + auto DME = Teuchos::rcp_dynamic_cast< + const Tempus_Test::DahlquistTestModel > (stepper->getModel()); + TEST_FLOATING_EQUALITY(DME->getLambda(), -1.0, relTol); + } + TEST_FLOATING_EQUALITY(dt, 1.0, relTol); + + const double x_0 = get_ele(*(x), 0); + const double xDot_0 = get_ele(*(xDot), 0); + TEST_FLOATING_EQUALITY(x_0, 1.0, relTol); // Should be x_0 + TEST_FLOATING_EQUALITY(xDot_0, -1.0, relTol); // Should be xDot_0 + TEST_ASSERT(std::abs(time) < relTol); + TEST_FLOATING_EQUALITY(dt, 1.0, relTol); + TEST_COMPARE(stageNumber, ==, -1); + break; + } + case StepperRKAppAction::BEGIN_STAGE: + case StepperRKAppAction::BEFORE_SOLVE: + { + const double X_i = get_ele(*(x), 0); + const double f_i = get_ele(*(xDot), 0); + + if (stageNumber == 0) { + TEST_FLOATING_EQUALITY(X_i, 1.0, relTol); + TEST_FLOATING_EQUALITY(time, 1.0, relTol); + TEST_ASSERT(std::abs(f_i) < relTol); + } else if (stageNumber == 1) { + TEST_FLOATING_EQUALITY(X_i, 0.5, relTol); + TEST_FLOATING_EQUALITY(f_i, -1.0, relTol); + TEST_ASSERT(std::abs(time) < relTol); + } else { + TEUCHOS_TEST_FOR_EXCEPT( !(-1 < stageNumber && stageNumber < 2)); + } + + break; + } + case StepperRKAppAction::AFTER_SOLVE: + case StepperRKAppAction::BEFORE_EXPLICIT_EVAL: + case StepperRKAppAction::END_STAGE: + { + const double X_i = get_ele(*(x), 0); // 0.5 + const double f_i = get_ele(*(xDot), 0); // -0.5 + + if (stageNumber == 0) { + // X_i = 1, f_1 = 0 + TEST_FLOATING_EQUALITY(X_i, 0.5, relTol); + TEST_FLOATING_EQUALITY(time, 1.0, relTol); + TEST_FLOATING_EQUALITY(f_i, -0.5, relTol); + } else if (stageNumber == 1) { + // X_i = , f_i = + TEST_FLOATING_EQUALITY(X_i, 0.75, relTol); + TEST_FLOATING_EQUALITY(f_i, -0.75, relTol); + TEST_ASSERT(std::abs(time) < relTol); + } else { + TEUCHOS_TEST_FOR_EXCEPT( !(-1 < stageNumber && stageNumber < 2)); + } + + break; + } + case StepperRKAppAction::END_STEP: + { + const double x_1 = get_ele(*(x), 0); + time = workingState->getTime(); + TEST_FLOATING_EQUALITY(x_1, 3.0/8.0, relTol); // Should be x_1 + TEST_FLOATING_EQUALITY(time, 1.0, relTol); + TEST_FLOATING_EQUALITY(dt, 1.0, relTol); + TEST_COMPARE(stageNumber, ==, -1); + + if (stepper->getUseEmbedded() == true) { + TEST_FLOATING_EQUALITY(workingState->getTolRel(), 1.0, relTol); + TEST_ASSERT(std::abs(workingState->getTolAbs()) < relTol); + // e = 0 from doxygen above. + TEST_ASSERT(std::abs(workingState->getErrorRel()) < relTol); + } + + } + + } + + + } + +private: + + Teuchos::FancyOStream & out; + bool & success; +}; + +// ************************************************************ +// ************************************************************ + +class StepperRKModifierEDIRK_TrapezoidaTest + : virtual public Tempus::StepperRKModifierBase +{ +public: + + /// Constructor + StepperRKModifierEDIRK_TrapezoidaTest(Teuchos::FancyOStream &Out, bool &Success) + : out(Out), success(Success) + {} + + // FSAL + /// Destructor + virtual ~StepperRKModifierEDIRK_TrapezoidaTest(){} + + /// Test the modify RK Stepper at the Action Locations. + virtual void modify( + Teuchos::RCP > sh, + Teuchos::RCP > stepper, + const typename Tempus::StepperRKAppAction::ACTION_LOCATION actLoc) + { + const double relTol = 1.0e-14; + auto stageNumber = stepper->getStageNumber(); + Teuchos::SerialDenseVector c = stepper->getTableau()->c(); + + auto currentState = sh->getCurrentState(); + auto workingState = sh->getWorkingState(); + const double dt = workingState->getTimeStep(); + double time = currentState->getTime(); + if (stageNumber >= 0) time += c(stageNumber)*dt; + + auto x = workingState->getX(); + auto xDot = workingState->getXDot(); + if (xDot == Teuchos::null) xDot = stepper->getStepperXDot(); + + switch (actLoc) { + case StepperRKAppAction::BEGIN_STEP: + { + { + auto DME = Teuchos::rcp_dynamic_cast< + const Tempus_Test::DahlquistTestModel > (stepper->getModel()); + TEST_FLOATING_EQUALITY(DME->getLambda(), -1.0, relTol); + } + TEST_FLOATING_EQUALITY(dt, 1.0, relTol); + + const double x_0 = get_ele(*(x), 0); + const double xDot_0 = get_ele(*(xDot), 0); + TEST_FLOATING_EQUALITY(x_0, 1.0, relTol); // Should be x_0 + TEST_FLOATING_EQUALITY(xDot_0, -1.0, relTol); // Should be xDot_0 + TEST_ASSERT(std::abs(time) < relTol); + TEST_FLOATING_EQUALITY(dt, 1.0, relTol); + TEST_COMPARE(stageNumber, ==, -1); + break; + } + case StepperRKAppAction::BEGIN_STAGE: + case StepperRKAppAction::BEFORE_SOLVE: + { + const double X_i = get_ele(*(x), 0); + const double f_i = get_ele(*(xDot), 0); + + if (stageNumber == 0) { + TEST_FLOATING_EQUALITY(X_i, 1.0, relTol); + TEST_ASSERT(std::abs(f_i) < relTol); + TEST_ASSERT(std::abs(time) < relTol); + } else if (stageNumber == 1) { + TEST_FLOATING_EQUALITY(X_i, 1.0, relTol); + TEST_FLOATING_EQUALITY(f_i, -1.0, relTol); + TEST_FLOATING_EQUALITY(time, 1.0, relTol); + } else { + TEUCHOS_TEST_FOR_EXCEPT( !(-1 < stageNumber && stageNumber < 2)); + } + + break; + } + case StepperRKAppAction::AFTER_SOLVE: + case StepperRKAppAction::BEFORE_EXPLICIT_EVAL: + case StepperRKAppAction::END_STAGE: + { + const double X_i = get_ele(*(x), 0); + const double f_i = get_ele(*(xDot), 0); + + if (stageNumber == 0) { + // X_i = 1, f_1 = 0 + TEST_FLOATING_EQUALITY(X_i, 1.0, relTol); + TEST_FLOATING_EQUALITY(f_i, -1.0, relTol); + TEST_ASSERT(std::abs(time) < relTol); + } else if (stageNumber == 1) { + // X_i = , f_i = + TEST_FLOATING_EQUALITY(X_i, 1.0/3.0, relTol); + TEST_FLOATING_EQUALITY(f_i, -1.0/3.0, relTol); + TEST_FLOATING_EQUALITY(time, 1.0, relTol); + } else { + TEUCHOS_TEST_FOR_EXCEPT( !(-1 < stageNumber && stageNumber < 2)); + } + + break; + } + case StepperRKAppAction::END_STEP: + { + const double x_1 = get_ele(*(x), 0); + time = workingState->getTime(); + TEST_FLOATING_EQUALITY(x_1, 1.0/3.0, relTol); // Should be x_1 + TEST_FLOATING_EQUALITY(time, 1.0, relTol); + TEST_FLOATING_EQUALITY(dt, 1.0, relTol); + TEST_COMPARE(stageNumber, ==, -1); + + if (stepper->getUseEmbedded() == true) { + TEST_FLOATING_EQUALITY(workingState->getTolRel(), 1.0, relTol); + TEST_ASSERT(std::abs(workingState->getTolAbs()) < relTol); + // e = 0 from doxygen above. + TEST_ASSERT(std::abs(workingState->getErrorRel()) < relTol); + } + + } + + } + + + } + +private: + + Teuchos::FancyOStream & out; + bool & success; +}; + +// ************************************************************ +// ************************************************************ +TEUCHOS_UNIT_TEST(DIRK, SDIRK21_Modifier) +{ + auto stepper = rcp(new Tempus::StepperSDIRK_21Pair()); + Teuchos::RCP > + model = rcp(new Tempus_Test::DahlquistTestModel()); + + auto modifier = rcp(new StepperRKModifierSDIRK21Test(out, success)); + + stepper->setModel(model); + stepper->setAppAction(modifier); + stepper->setICConsistency("Consistent"); + stepper->setUseFSAL(false); + stepper->initialize(); + + // Create a SolutionHistory. + auto solutionHistory = Tempus::createSolutionHistoryME(model); + + // Take one time step. + stepper->setInitialConditions(solutionHistory); + solutionHistory->initWorkingState(); + double dt = 1.0; + solutionHistory->getWorkingState()->setTimeStep(dt); + solutionHistory->getWorkingState()->setTime(dt); + stepper->takeStep(solutionHistory); // Primary testing occurs in + // modifierX during takeStep(). + // Test stepper properties. + TEUCHOS_ASSERT(stepper->getOrder() == 2); +} + + +// ************************************************************ +// ************************************************************ +TEUCHOS_UNIT_TEST(DIRK, EDIRK_Trapezoida_Modifier) +{ + auto stepper = rcp(new Tempus::StepperEDIRK_TrapezoidalRule()); + Teuchos::RCP > + model = rcp(new Tempus_Test::DahlquistTestModel()); + + auto modifier = rcp(new StepperRKModifierEDIRK_TrapezoidaTest(out, success)); + + stepper->setModel(model); + stepper->setAppAction(modifier); + stepper->setICConsistency("Consistent"); + stepper->setUseFSAL(false); + stepper->initialize(); + + // Create a SolutionHistory. + auto solutionHistory = Tempus::createSolutionHistoryME(model); + + // Take one time step. + stepper->setInitialConditions(solutionHistory); + solutionHistory->initWorkingState(); + double dt = 1.0; + solutionHistory->getWorkingState()->setTimeStep(dt); + solutionHistory->getWorkingState()->setTime(dt); + stepper->takeStep(solutionHistory); // Primary testing occurs in + // modifierX during takeStep(). + // Test stepper properties. + TEUCHOS_ASSERT(stepper->getOrder() == 2); +} + +} // namespace Tempus_Test From 29f6d76f88197e9e90dc46e57314263d1879a69b Mon Sep 17 00:00:00 2001 From: Sidafa Conde Date: Sat, 27 Feb 2021 09:53:13 -0500 Subject: [PATCH 2/2] Tempus: split UnitTest_DIRK - move SDIRK21_Modifier to Tempus_UnitTest_SDIRK_21Pair - move EDIRK_Trapezoida_Modifier to Tempus_UnitTest_EDIRK_TrapezoidalRule - and remove the unnecessary Tempus_UnitTest_DIRK file since it is no longer needed --- packages/tempus/unit_test/CMakeLists.txt | 7 - .../tempus/unit_test/Tempus_UnitTest_DIRK.cpp | 343 ------------------ .../Tempus_UnitTest_EDIRK_TrapezoidalRule.cpp | 162 ++++++++- .../Tempus_UnitTest_SDIRK_21Pair.cpp | 161 +++++++- 4 files changed, 321 insertions(+), 352 deletions(-) delete mode 100644 packages/tempus/unit_test/Tempus_UnitTest_DIRK.cpp diff --git a/packages/tempus/unit_test/CMakeLists.txt b/packages/tempus/unit_test/CMakeLists.txt index 931aa47ed89a..8d79cca20684 100644 --- a/packages/tempus/unit_test/CMakeLists.txt +++ b/packages/tempus/unit_test/CMakeLists.txt @@ -414,13 +414,6 @@ TRIBITS_ADD_EXECUTABLE_AND_TEST( NUM_MPI_PROCS 1 ) -TRIBITS_ADD_EXECUTABLE_AND_TEST( - UnitTest_DIRK - SOURCES Tempus_UnitTest_DIRK.cpp ${TEUCHOS_STD_UNIT_TEST_MAIN} - TESTONLYLIBS tempus_test_models - NUM_MPI_PROCS 1 - ) - TRIBITS_COPY_FILES_TO_BINARY_DIR(UnitTest_NewmarkImplicitAForm_CopyFiles DEST_FILES Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder.xml Tempus_NewmarkExplicitAForm_HarmonicOscillator_Damped.xml EXEDEPS UnitTest_NewmarkImplicitAForm diff --git a/packages/tempus/unit_test/Tempus_UnitTest_DIRK.cpp b/packages/tempus/unit_test/Tempus_UnitTest_DIRK.cpp deleted file mode 100644 index 352cc55a94e3..000000000000 --- a/packages/tempus/unit_test/Tempus_UnitTest_DIRK.cpp +++ /dev/null @@ -1,343 +0,0 @@ -// @HEADER -// **************************************************************************** -// Tempus: Copyright (2017) Sandia Corporation -// -// Distributed under BSD 3-clause license (See accompanying file Copyright.txt) -// **************************************************************************** -// @HEADER - -#include "Teuchos_UnitTestHarness.hpp" - -#include "../TestModels/DahlquistTestModel.hpp" - -#include "Tempus_SolutionHistory.hpp" -#include "Tempus_StepperRKButcherTableau.hpp" -#include "Tempus_StepperRKBase.hpp" -#include "Tempus_StepperRKModifierBase.hpp" -#include "Tempus_StepperRKModifierXBase.hpp" - - -namespace Tempus_Unit_Test { - -using Teuchos::RCP; -using Teuchos::rcp; -using Teuchos::rcp_const_cast; -using Teuchos::rcp_dynamic_cast; - - -class StepperRKModifierSDIRK21Test - : virtual public Tempus::StepperRKModifierBase -{ -public: - - /// Constructor - StepperRKModifierSDIRK21Test(Teuchos::FancyOStream &Out, bool &Success) - : out(Out), success(Success) - {} - - /// Destructor - virtual ~StepperRKModifierSDIRK21Test(){} - - /// Test the modify RK Stepper at the Action Locations. - virtual void modify( - Teuchos::RCP > sh, - Teuchos::RCP > stepper, - const typename Tempus::StepperRKAppAction::ACTION_LOCATION actLoc) - { - const double relTol = 1.0e-14; - auto stageNumber = stepper->getStageNumber(); - Teuchos::SerialDenseVector c = stepper->getTableau()->c(); - - auto currentState = sh->getCurrentState(); - auto workingState = sh->getWorkingState(); - const double dt = workingState->getTimeStep(); - double time = currentState->getTime(); - if (stageNumber >= 0) time += c(stageNumber)*dt; - - auto x = workingState->getX(); - auto xDot = workingState->getXDot(); - if (xDot == Teuchos::null) xDot = stepper->getStepperXDot(); - - switch (actLoc) { - case StepperRKAppAction::BEGIN_STEP: - { - { - auto DME = Teuchos::rcp_dynamic_cast< - const Tempus_Test::DahlquistTestModel > (stepper->getModel()); - TEST_FLOATING_EQUALITY(DME->getLambda(), -1.0, relTol); - } - TEST_FLOATING_EQUALITY(dt, 1.0, relTol); - - const double x_0 = get_ele(*(x), 0); - const double xDot_0 = get_ele(*(xDot), 0); - TEST_FLOATING_EQUALITY(x_0, 1.0, relTol); // Should be x_0 - TEST_FLOATING_EQUALITY(xDot_0, -1.0, relTol); // Should be xDot_0 - TEST_ASSERT(std::abs(time) < relTol); - TEST_FLOATING_EQUALITY(dt, 1.0, relTol); - TEST_COMPARE(stageNumber, ==, -1); - break; - } - case StepperRKAppAction::BEGIN_STAGE: - case StepperRKAppAction::BEFORE_SOLVE: - { - const double X_i = get_ele(*(x), 0); - const double f_i = get_ele(*(xDot), 0); - - if (stageNumber == 0) { - TEST_FLOATING_EQUALITY(X_i, 1.0, relTol); - TEST_FLOATING_EQUALITY(time, 1.0, relTol); - TEST_ASSERT(std::abs(f_i) < relTol); - } else if (stageNumber == 1) { - TEST_FLOATING_EQUALITY(X_i, 0.5, relTol); - TEST_FLOATING_EQUALITY(f_i, -1.0, relTol); - TEST_ASSERT(std::abs(time) < relTol); - } else { - TEUCHOS_TEST_FOR_EXCEPT( !(-1 < stageNumber && stageNumber < 2)); - } - - break; - } - case StepperRKAppAction::AFTER_SOLVE: - case StepperRKAppAction::BEFORE_EXPLICIT_EVAL: - case StepperRKAppAction::END_STAGE: - { - const double X_i = get_ele(*(x), 0); // 0.5 - const double f_i = get_ele(*(xDot), 0); // -0.5 - - if (stageNumber == 0) { - // X_i = 1, f_1 = 0 - TEST_FLOATING_EQUALITY(X_i, 0.5, relTol); - TEST_FLOATING_EQUALITY(time, 1.0, relTol); - TEST_FLOATING_EQUALITY(f_i, -0.5, relTol); - } else if (stageNumber == 1) { - // X_i = , f_i = - TEST_FLOATING_EQUALITY(X_i, 0.75, relTol); - TEST_FLOATING_EQUALITY(f_i, -0.75, relTol); - TEST_ASSERT(std::abs(time) < relTol); - } else { - TEUCHOS_TEST_FOR_EXCEPT( !(-1 < stageNumber && stageNumber < 2)); - } - - break; - } - case StepperRKAppAction::END_STEP: - { - const double x_1 = get_ele(*(x), 0); - time = workingState->getTime(); - TEST_FLOATING_EQUALITY(x_1, 3.0/8.0, relTol); // Should be x_1 - TEST_FLOATING_EQUALITY(time, 1.0, relTol); - TEST_FLOATING_EQUALITY(dt, 1.0, relTol); - TEST_COMPARE(stageNumber, ==, -1); - - if (stepper->getUseEmbedded() == true) { - TEST_FLOATING_EQUALITY(workingState->getTolRel(), 1.0, relTol); - TEST_ASSERT(std::abs(workingState->getTolAbs()) < relTol); - // e = 0 from doxygen above. - TEST_ASSERT(std::abs(workingState->getErrorRel()) < relTol); - } - - } - - } - - - } - -private: - - Teuchos::FancyOStream & out; - bool & success; -}; - -// ************************************************************ -// ************************************************************ - -class StepperRKModifierEDIRK_TrapezoidaTest - : virtual public Tempus::StepperRKModifierBase -{ -public: - - /// Constructor - StepperRKModifierEDIRK_TrapezoidaTest(Teuchos::FancyOStream &Out, bool &Success) - : out(Out), success(Success) - {} - - // FSAL - /// Destructor - virtual ~StepperRKModifierEDIRK_TrapezoidaTest(){} - - /// Test the modify RK Stepper at the Action Locations. - virtual void modify( - Teuchos::RCP > sh, - Teuchos::RCP > stepper, - const typename Tempus::StepperRKAppAction::ACTION_LOCATION actLoc) - { - const double relTol = 1.0e-14; - auto stageNumber = stepper->getStageNumber(); - Teuchos::SerialDenseVector c = stepper->getTableau()->c(); - - auto currentState = sh->getCurrentState(); - auto workingState = sh->getWorkingState(); - const double dt = workingState->getTimeStep(); - double time = currentState->getTime(); - if (stageNumber >= 0) time += c(stageNumber)*dt; - - auto x = workingState->getX(); - auto xDot = workingState->getXDot(); - if (xDot == Teuchos::null) xDot = stepper->getStepperXDot(); - - switch (actLoc) { - case StepperRKAppAction::BEGIN_STEP: - { - { - auto DME = Teuchos::rcp_dynamic_cast< - const Tempus_Test::DahlquistTestModel > (stepper->getModel()); - TEST_FLOATING_EQUALITY(DME->getLambda(), -1.0, relTol); - } - TEST_FLOATING_EQUALITY(dt, 1.0, relTol); - - const double x_0 = get_ele(*(x), 0); - const double xDot_0 = get_ele(*(xDot), 0); - TEST_FLOATING_EQUALITY(x_0, 1.0, relTol); // Should be x_0 - TEST_FLOATING_EQUALITY(xDot_0, -1.0, relTol); // Should be xDot_0 - TEST_ASSERT(std::abs(time) < relTol); - TEST_FLOATING_EQUALITY(dt, 1.0, relTol); - TEST_COMPARE(stageNumber, ==, -1); - break; - } - case StepperRKAppAction::BEGIN_STAGE: - case StepperRKAppAction::BEFORE_SOLVE: - { - const double X_i = get_ele(*(x), 0); - const double f_i = get_ele(*(xDot), 0); - - if (stageNumber == 0) { - TEST_FLOATING_EQUALITY(X_i, 1.0, relTol); - TEST_ASSERT(std::abs(f_i) < relTol); - TEST_ASSERT(std::abs(time) < relTol); - } else if (stageNumber == 1) { - TEST_FLOATING_EQUALITY(X_i, 1.0, relTol); - TEST_FLOATING_EQUALITY(f_i, -1.0, relTol); - TEST_FLOATING_EQUALITY(time, 1.0, relTol); - } else { - TEUCHOS_TEST_FOR_EXCEPT( !(-1 < stageNumber && stageNumber < 2)); - } - - break; - } - case StepperRKAppAction::AFTER_SOLVE: - case StepperRKAppAction::BEFORE_EXPLICIT_EVAL: - case StepperRKAppAction::END_STAGE: - { - const double X_i = get_ele(*(x), 0); - const double f_i = get_ele(*(xDot), 0); - - if (stageNumber == 0) { - // X_i = 1, f_1 = 0 - TEST_FLOATING_EQUALITY(X_i, 1.0, relTol); - TEST_FLOATING_EQUALITY(f_i, -1.0, relTol); - TEST_ASSERT(std::abs(time) < relTol); - } else if (stageNumber == 1) { - // X_i = , f_i = - TEST_FLOATING_EQUALITY(X_i, 1.0/3.0, relTol); - TEST_FLOATING_EQUALITY(f_i, -1.0/3.0, relTol); - TEST_FLOATING_EQUALITY(time, 1.0, relTol); - } else { - TEUCHOS_TEST_FOR_EXCEPT( !(-1 < stageNumber && stageNumber < 2)); - } - - break; - } - case StepperRKAppAction::END_STEP: - { - const double x_1 = get_ele(*(x), 0); - time = workingState->getTime(); - TEST_FLOATING_EQUALITY(x_1, 1.0/3.0, relTol); // Should be x_1 - TEST_FLOATING_EQUALITY(time, 1.0, relTol); - TEST_FLOATING_EQUALITY(dt, 1.0, relTol); - TEST_COMPARE(stageNumber, ==, -1); - - if (stepper->getUseEmbedded() == true) { - TEST_FLOATING_EQUALITY(workingState->getTolRel(), 1.0, relTol); - TEST_ASSERT(std::abs(workingState->getTolAbs()) < relTol); - // e = 0 from doxygen above. - TEST_ASSERT(std::abs(workingState->getErrorRel()) < relTol); - } - - } - - } - - - } - -private: - - Teuchos::FancyOStream & out; - bool & success; -}; - -// ************************************************************ -// ************************************************************ -TEUCHOS_UNIT_TEST(DIRK, SDIRK21_Modifier) -{ - auto stepper = rcp(new Tempus::StepperSDIRK_21Pair()); - Teuchos::RCP > - model = rcp(new Tempus_Test::DahlquistTestModel()); - - auto modifier = rcp(new StepperRKModifierSDIRK21Test(out, success)); - - stepper->setModel(model); - stepper->setAppAction(modifier); - stepper->setICConsistency("Consistent"); - stepper->setUseFSAL(false); - stepper->initialize(); - - // Create a SolutionHistory. - auto solutionHistory = Tempus::createSolutionHistoryME(model); - - // Take one time step. - stepper->setInitialConditions(solutionHistory); - solutionHistory->initWorkingState(); - double dt = 1.0; - solutionHistory->getWorkingState()->setTimeStep(dt); - solutionHistory->getWorkingState()->setTime(dt); - stepper->takeStep(solutionHistory); // Primary testing occurs in - // modifierX during takeStep(). - // Test stepper properties. - TEUCHOS_ASSERT(stepper->getOrder() == 2); -} - - -// ************************************************************ -// ************************************************************ -TEUCHOS_UNIT_TEST(DIRK, EDIRK_Trapezoida_Modifier) -{ - auto stepper = rcp(new Tempus::StepperEDIRK_TrapezoidalRule()); - Teuchos::RCP > - model = rcp(new Tempus_Test::DahlquistTestModel()); - - auto modifier = rcp(new StepperRKModifierEDIRK_TrapezoidaTest(out, success)); - - stepper->setModel(model); - stepper->setAppAction(modifier); - stepper->setICConsistency("Consistent"); - stepper->setUseFSAL(false); - stepper->initialize(); - - // Create a SolutionHistory. - auto solutionHistory = Tempus::createSolutionHistoryME(model); - - // Take one time step. - stepper->setInitialConditions(solutionHistory); - solutionHistory->initWorkingState(); - double dt = 1.0; - solutionHistory->getWorkingState()->setTimeStep(dt); - solutionHistory->getWorkingState()->setTime(dt); - stepper->takeStep(solutionHistory); // Primary testing occurs in - // modifierX during takeStep(). - // Test stepper properties. - TEUCHOS_ASSERT(stepper->getOrder() == 2); -} - -} // namespace Tempus_Test diff --git a/packages/tempus/unit_test/Tempus_UnitTest_EDIRK_TrapezoidalRule.cpp b/packages/tempus/unit_test/Tempus_UnitTest_EDIRK_TrapezoidalRule.cpp index 571b0c3c3b3f..854116e75bc0 100644 --- a/packages/tempus/unit_test/Tempus_UnitTest_EDIRK_TrapezoidalRule.cpp +++ b/packages/tempus/unit_test/Tempus_UnitTest_EDIRK_TrapezoidalRule.cpp @@ -16,9 +16,10 @@ #include "Tempus_UnitTest_Utils.hpp" #include "../TestModels/SinCosModel.hpp" -#include "../TestModels/VanDerPolModel.hpp" +#include "../TestModels/DahlquistTestModel.hpp" #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp" +#include "Tempus_StepperRKModifierBase.hpp" #include #include @@ -64,4 +65,163 @@ TEUCHOS_UNIT_TEST(EDIRK_TrapezoidalRule, AppAction) } +// ************************************************************ +// ************************************************************ + +class StepperRKModifierEDIRK_TrapezoidaTest + : virtual public Tempus::StepperRKModifierBase +{ +public: + + /// Constructor + StepperRKModifierEDIRK_TrapezoidaTest(Teuchos::FancyOStream &Out, bool &Success) + : out(Out), success(Success) + {} + + // FSAL + /// Destructor + virtual ~StepperRKModifierEDIRK_TrapezoidaTest(){} + + /// Test the modify RK Stepper at the Action Locations. + virtual void modify( + Teuchos::RCP > sh, + Teuchos::RCP > stepper, + const typename Tempus::StepperRKAppAction::ACTION_LOCATION actLoc) + { + const double relTol = 1.0e-14; + auto stageNumber = stepper->getStageNumber(); + Teuchos::SerialDenseVector c = stepper->getTableau()->c(); + + auto currentState = sh->getCurrentState(); + auto workingState = sh->getWorkingState(); + const double dt = workingState->getTimeStep(); + double time = currentState->getTime(); + if (stageNumber >= 0) time += c(stageNumber)*dt; + + auto x = workingState->getX(); + auto xDot = workingState->getXDot(); + if (xDot == Teuchos::null) xDot = stepper->getStepperXDot(); + + switch (actLoc) { + case StepperRKAppAction::BEGIN_STEP: + { + { + auto DME = Teuchos::rcp_dynamic_cast< + const Tempus_Test::DahlquistTestModel > (stepper->getModel()); + TEST_FLOATING_EQUALITY(DME->getLambda(), -1.0, relTol); + } + TEST_FLOATING_EQUALITY(dt, 1.0, relTol); + + const double x_0 = get_ele(*(x), 0); + const double xDot_0 = get_ele(*(xDot), 0); + TEST_FLOATING_EQUALITY(x_0, 1.0, relTol); // Should be x_0 + TEST_FLOATING_EQUALITY(xDot_0, -1.0, relTol); // Should be xDot_0 + TEST_ASSERT(std::abs(time) < relTol); + TEST_FLOATING_EQUALITY(dt, 1.0, relTol); + TEST_COMPARE(stageNumber, ==, -1); + break; + } + case StepperRKAppAction::BEGIN_STAGE: + case StepperRKAppAction::BEFORE_SOLVE: + { + const double X_i = get_ele(*(x), 0); + const double f_i = get_ele(*(xDot), 0); + + if (stageNumber == 0) { + TEST_FLOATING_EQUALITY(X_i, 1.0, relTol); + TEST_ASSERT(std::abs(f_i) < relTol); + TEST_ASSERT(std::abs(time) < relTol); + } else if (stageNumber == 1) { + TEST_FLOATING_EQUALITY(X_i, 1.0, relTol); + TEST_FLOATING_EQUALITY(f_i, -1.0, relTol); + TEST_FLOATING_EQUALITY(time, 1.0, relTol); + } else { + TEUCHOS_TEST_FOR_EXCEPT( !(-1 < stageNumber && stageNumber < 2)); + } + + break; + } + case StepperRKAppAction::AFTER_SOLVE: + case StepperRKAppAction::BEFORE_EXPLICIT_EVAL: + case StepperRKAppAction::END_STAGE: + { + const double X_i = get_ele(*(x), 0); + const double f_i = get_ele(*(xDot), 0); + + if (stageNumber == 0) { + // X_i = 1, f_1 = 0 + TEST_FLOATING_EQUALITY(X_i, 1.0, relTol); + TEST_FLOATING_EQUALITY(f_i, -1.0, relTol); + TEST_ASSERT(std::abs(time) < relTol); + } else if (stageNumber == 1) { + // X_i = , f_i = + TEST_FLOATING_EQUALITY(X_i, 1.0/3.0, relTol); + TEST_FLOATING_EQUALITY(f_i, -1.0/3.0, relTol); + TEST_FLOATING_EQUALITY(time, 1.0, relTol); + } else { + TEUCHOS_TEST_FOR_EXCEPT( !(-1 < stageNumber && stageNumber < 2)); + } + + break; + } + case StepperRKAppAction::END_STEP: + { + const double x_1 = get_ele(*(x), 0); + time = workingState->getTime(); + TEST_FLOATING_EQUALITY(x_1, 1.0/3.0, relTol); // Should be x_1 + TEST_FLOATING_EQUALITY(time, 1.0, relTol); + TEST_FLOATING_EQUALITY(dt, 1.0, relTol); + TEST_COMPARE(stageNumber, ==, -1); + + if (stepper->getUseEmbedded() == true) { + TEST_FLOATING_EQUALITY(workingState->getTolRel(), 1.0, relTol); + TEST_ASSERT(std::abs(workingState->getTolAbs()) < relTol); + // e = 0 from doxygen above. + TEST_ASSERT(std::abs(workingState->getErrorRel()) < relTol); + } + + } + + } + + + } + +private: + + Teuchos::FancyOStream & out; + bool & success; +}; + + +// ************************************************************ +// ************************************************************ +TEUCHOS_UNIT_TEST(DIRK, EDIRK_Trapezoida_Modifier) +{ + auto stepper = rcp(new Tempus::StepperEDIRK_TrapezoidalRule()); + Teuchos::RCP > + model = rcp(new Tempus_Test::DahlquistTestModel()); + + auto modifier = rcp(new StepperRKModifierEDIRK_TrapezoidaTest(out, success)); + + stepper->setModel(model); + stepper->setAppAction(modifier); + stepper->setICConsistency("Consistent"); + stepper->setUseFSAL(false); + stepper->initialize(); + + // Create a SolutionHistory. + auto solutionHistory = Tempus::createSolutionHistoryME(model); + + // Take one time step. + stepper->setInitialConditions(solutionHistory); + solutionHistory->initWorkingState(); + double dt = 1.0; + solutionHistory->getWorkingState()->setTimeStep(dt); + solutionHistory->getWorkingState()->setTime(dt); + stepper->takeStep(solutionHistory); // Primary testing occurs in + // modifierX during takeStep(). + // Test stepper properties. + TEUCHOS_ASSERT(stepper->getOrder() == 2); +} } // namespace Tempus_Test diff --git a/packages/tempus/unit_test/Tempus_UnitTest_SDIRK_21Pair.cpp b/packages/tempus/unit_test/Tempus_UnitTest_SDIRK_21Pair.cpp index ba5d4d970370..eaf9f0c05230 100644 --- a/packages/tempus/unit_test/Tempus_UnitTest_SDIRK_21Pair.cpp +++ b/packages/tempus/unit_test/Tempus_UnitTest_SDIRK_21Pair.cpp @@ -16,9 +16,10 @@ #include "Tempus_UnitTest_Utils.hpp" #include "../TestModels/SinCosModel.hpp" -#include "../TestModels/VanDerPolModel.hpp" +#include "../TestModels/DahlquistTestModel.hpp" #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp" +#include "Tempus_StepperRKModifierBase.hpp" #include #include @@ -64,4 +65,162 @@ TEUCHOS_UNIT_TEST(SDIRK_21Pair, AppAction) } +// ************************************************************ +// ************************************************************ +class StepperRKModifierSDIRK21Test + : virtual public Tempus::StepperRKModifierBase +{ +public: + + /// Constructor + StepperRKModifierSDIRK21Test(Teuchos::FancyOStream &Out, bool &Success) + : out(Out), success(Success) + {} + + /// Destructor + virtual ~StepperRKModifierSDIRK21Test(){} + + /// Test the modify RK Stepper at the Action Locations. + virtual void modify( + Teuchos::RCP > sh, + Teuchos::RCP > stepper, + const typename Tempus::StepperRKAppAction::ACTION_LOCATION actLoc) + { + const double relTol = 1.0e-14; + auto stageNumber = stepper->getStageNumber(); + Teuchos::SerialDenseVector c = stepper->getTableau()->c(); + + auto currentState = sh->getCurrentState(); + auto workingState = sh->getWorkingState(); + const double dt = workingState->getTimeStep(); + double time = currentState->getTime(); + if (stageNumber >= 0) time += c(stageNumber)*dt; + + auto x = workingState->getX(); + auto xDot = workingState->getXDot(); + if (xDot == Teuchos::null) xDot = stepper->getStepperXDot(); + + switch (actLoc) { + case StepperRKAppAction::BEGIN_STEP: + { + { + auto DME = Teuchos::rcp_dynamic_cast< + const Tempus_Test::DahlquistTestModel > (stepper->getModel()); + TEST_FLOATING_EQUALITY(DME->getLambda(), -1.0, relTol); + } + TEST_FLOATING_EQUALITY(dt, 1.0, relTol); + + const double x_0 = get_ele(*(x), 0); + const double xDot_0 = get_ele(*(xDot), 0); + TEST_FLOATING_EQUALITY(x_0, 1.0, relTol); // Should be x_0 + TEST_FLOATING_EQUALITY(xDot_0, -1.0, relTol); // Should be xDot_0 + TEST_ASSERT(std::abs(time) < relTol); + TEST_FLOATING_EQUALITY(dt, 1.0, relTol); + TEST_COMPARE(stageNumber, ==, -1); + break; + } + case StepperRKAppAction::BEGIN_STAGE: + case StepperRKAppAction::BEFORE_SOLVE: + { + const double X_i = get_ele(*(x), 0); + const double f_i = get_ele(*(xDot), 0); + + if (stageNumber == 0) { + TEST_FLOATING_EQUALITY(X_i, 1.0, relTol); + TEST_FLOATING_EQUALITY(time, 1.0, relTol); + TEST_ASSERT(std::abs(f_i) < relTol); + } else if (stageNumber == 1) { + TEST_FLOATING_EQUALITY(X_i, 0.5, relTol); + TEST_FLOATING_EQUALITY(f_i, -1.0, relTol); + TEST_ASSERT(std::abs(time) < relTol); + } else { + TEUCHOS_TEST_FOR_EXCEPT( !(-1 < stageNumber && stageNumber < 2)); + } + + break; + } + case StepperRKAppAction::AFTER_SOLVE: + case StepperRKAppAction::BEFORE_EXPLICIT_EVAL: + case StepperRKAppAction::END_STAGE: + { + const double X_i = get_ele(*(x), 0); // 0.5 + const double f_i = get_ele(*(xDot), 0); // -0.5 + + if (stageNumber == 0) { + // X_i = 1, f_1 = 0 + TEST_FLOATING_EQUALITY(X_i, 0.5, relTol); + TEST_FLOATING_EQUALITY(time, 1.0, relTol); + TEST_FLOATING_EQUALITY(f_i, -0.5, relTol); + } else if (stageNumber == 1) { + // X_i = , f_i = + TEST_FLOATING_EQUALITY(X_i, 0.75, relTol); + TEST_FLOATING_EQUALITY(f_i, -0.75, relTol); + TEST_ASSERT(std::abs(time) < relTol); + } else { + TEUCHOS_TEST_FOR_EXCEPT( !(-1 < stageNumber && stageNumber < 2)); + } + + break; + } + case StepperRKAppAction::END_STEP: + { + const double x_1 = get_ele(*(x), 0); + time = workingState->getTime(); + TEST_FLOATING_EQUALITY(x_1, 3.0/8.0, relTol); // Should be x_1 + TEST_FLOATING_EQUALITY(time, 1.0, relTol); + TEST_FLOATING_EQUALITY(dt, 1.0, relTol); + TEST_COMPARE(stageNumber, ==, -1); + + if (stepper->getUseEmbedded() == true) { + TEST_FLOATING_EQUALITY(workingState->getTolRel(), 1.0, relTol); + TEST_ASSERT(std::abs(workingState->getTolAbs()) < relTol); + // e = 0 from doxygen above. + TEST_ASSERT(std::abs(workingState->getErrorRel()) < relTol); + } + + } + + } + + + } + +private: + + Teuchos::FancyOStream & out; + bool & success; +}; + + +// ************************************************************ +// ************************************************************ +TEUCHOS_UNIT_TEST(SDIRK_21Pair, Modifier) +{ + auto stepper = rcp(new Tempus::StepperSDIRK_21Pair()); + Teuchos::RCP > + model = rcp(new Tempus_Test::DahlquistTestModel()); + + auto modifier = rcp(new StepperRKModifierSDIRK21Test(out, success)); + + stepper->setModel(model); + stepper->setAppAction(modifier); + stepper->setICConsistency("Consistent"); + stepper->setUseFSAL(false); + stepper->initialize(); + + // Create a SolutionHistory. + auto solutionHistory = Tempus::createSolutionHistoryME(model); + + // Take one time step. + stepper->setInitialConditions(solutionHistory); + solutionHistory->initWorkingState(); + double dt = 1.0; + solutionHistory->getWorkingState()->setTimeStep(dt); + solutionHistory->getWorkingState()->setTime(dt); + stepper->takeStep(solutionHistory); // Primary testing occurs in + // modifierX during takeStep(). + // Test stepper properties. + TEUCHOS_ASSERT(stepper->getOrder() == 2); +} + } // namespace Tempus_Test