From a81d9d866a03cda003fa1d146586c0f362bf036a Mon Sep 17 00:00:00 2001 From: Jacob Domagala Date: Mon, 25 Sep 2023 11:18:58 +0200 Subject: [PATCH] Tempus: Use CDR_Model_Tpetra for BDF2 Test --- packages/tempus/test/BDF2/Tempus_BDF2Test.cpp | 107 ++++++++++-------- 1 file changed, 62 insertions(+), 45 deletions(-) diff --git a/packages/tempus/test/BDF2/Tempus_BDF2Test.cpp b/packages/tempus/test/BDF2/Tempus_BDF2Test.cpp index 06fa5b3af6e8..cdfab2b5ea28 100644 --- a/packages/tempus/test/BDF2/Tempus_BDF2Test.cpp +++ b/packages/tempus/test/BDF2/Tempus_BDF2Test.cpp @@ -23,6 +23,7 @@ #include "Stratimikos_DefaultLinearSolverBuilder.hpp" #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp" +#include "Tpetra_Core.hpp" #ifdef Tempus_ENABLE_MPI #include "Epetra_MpiComm.h" @@ -475,46 +476,34 @@ TEUCHOS_UNIT_TEST(BDF2, SinCosAdapt) } -// ************************************************************ -// ************************************************************ -TEUCHOS_UNIT_TEST(BDF2, CDR) -{ - // Create a communicator for Epetra objects - RCP comm; -#ifdef Tempus_ENABLE_MPI - comm = rcp(new Epetra_MpiComm(MPI_COMM_WORLD)); -#else - comm = rcp(new Epetra_SerialComm); -#endif - +template +void CDR_Test(const Comm& comm, const int commSize, Teuchos::FancyOStream &out, bool& success){ RCP > integrator; std::vector>> solutions; std::vector>> solutionsDot; std::vector StepSize; // Read params from .xml file - RCP pList = - getParametersFromXmlFile("Tempus_BDF2_CDR.xml"); + auto pList = getParametersFromXmlFile("Tempus_BDF2_CDR.xml"); //Set initial time step = 2*dt specified in input file (for convergence study) // - RCP pl = sublist(pList, "Tempus", true); - double dt = pl->sublist("Demo Integrator") - .sublist("Time Step Control").get("Initial Time Step"); + auto pl = sublist(pList, "Tempus", true); + auto dt = pl->sublist("Demo Integrator").sublist("Time Step Control").get("Initial Time Step"); dt *= 2.0; - RCP model_pl = sublist(pList, "CDR Model", true); + auto model_pl = sublist(pList, "CDR Model", true); - const int nTimeStepSizes = model_pl->get("Number of Time Step Sizes", 5); + const auto nTimeStepSizes = model_pl->get("Number of Time Step Sizes", 5); for (int n=0; nget("num elements"); - const double left_end = model_pl->get("left end"); - const double right_end = model_pl->get("right end"); - const double a_convection = model_pl->get("a (convection)"); - const double k_source = model_pl->get("k (source)"); + const auto left_end = model_pl->get("left end"); + const auto right_end = model_pl->get("right end"); + const auto a_convection = model_pl->get("a (convection)"); + const auto k_source = model_pl->get("k (source)"); - auto model = rcp(new Tempus_Test::CDR_Model(comm, + auto model = rcp(new Model(comm, num_elements, left_end, right_end, @@ -529,8 +518,7 @@ TEUCHOS_UNIT_TEST(BDF2, CDR) p->set("Preconditioner Type", "None"); builder.setParameterList(p); - RCP< ::Thyra::LinearOpWithSolveFactoryBase > - lowsFactory = builder.createLinearSolveStrategy(""); + auto lowsFactory = builder.createLinearSolveStrategy(""); model->set_W_factory(lowsFactory); @@ -564,23 +552,23 @@ TEUCHOS_UNIT_TEST(BDF2, CDR) // Output finest temporal solution for plotting // This only works for ONE MPI process - if ((n == nTimeStepSizes-1) && (comm->NumProc() == 1)) { + if ((n == nTimeStepSizes-1) && (commSize == 1)) { std::ofstream ftmp("Tempus_BDF2_CDR.dat"); ftmp << "TITLE=\"BDF2 Solution to CDR\"\n" << "VARIABLES=\"z\",\"T\"\n"; - const double dx = std::fabs(left_end-right_end) / + const auto dx = std::fabs(left_end-right_end) / static_cast(num_elements); - RCP > solutionHistory = + auto solutionHistory = integrator->getSolutionHistory(); int nStates = solutionHistory->getNumStates(); for (int i=0; i > solutionState = (*solutionHistory)[i]; - RCP > x = solutionState->getX(); - double ttime = solutionState->getTime(); + auto solutionState = (*solutionHistory)[i]; + auto x = solutionState->getX(); + auto ttime = solutionState->getTime(); ftmp << "ZONE T=\"Time="<(j) * dx; + const auto x_coord = left_end + static_cast(j) * dx; ftmp << x_coord << " "; } ftmp << std::endl; @@ -597,8 +585,8 @@ TEUCHOS_UNIT_TEST(BDF2, CDR) double xDotSlope = 0.0; std::vector xErrorNorm; std::vector xDotErrorNorm; - RCP > stepper = integrator->getStepper(); - double order = stepper->getOrder(); + auto stepper = integrator->getStepper(); + auto order = stepper->getOrder(); writeOrderError("Tempus_BDF2_CDR-Error.dat", stepper, StepSize, solutions, xErrorNorm, xSlope, @@ -614,21 +602,20 @@ TEUCHOS_UNIT_TEST(BDF2, CDR) // Write fine mesh solution at final time // This only works for ONE MPI process - if (comm->NumProc() == 1) { - RCP pListCDR = - getParametersFromXmlFile("Tempus_BDF2_CDR.xml"); - RCP model_pl_CDR = sublist(pListCDR, "CDR Model", true); - const int num_elements = model_pl_CDR->get("num elements"); - const double left_end = model_pl_CDR->get("left end"); - const double right_end = model_pl_CDR->get("right end"); + if (commSize == 1) { + auto pListCDR = getParametersFromXmlFile("Tempus_BDF2_CDR.xml"); + auto model_pl_CDR = sublist(pListCDR, "CDR Model", true); + const auto num_elements = model_pl_CDR->get("num elements"); + const auto left_end = model_pl_CDR->get("left end"); + const auto right_end = model_pl_CDR->get("right end"); - const Thyra::VectorBase& x = *(solutions[solutions.size()-1]); + const auto& x = *(solutions[solutions.size()-1]); std::ofstream ftmp("Tempus_BDF2_CDR-Solution.dat"); for (int n = 0; n < num_elements+1; n++) { - const double dx = std::fabs(left_end-right_end) / + const auto dx = std::fabs(left_end-right_end) / static_cast(num_elements); - const double x_coord = left_end + static_cast(n) * dx; + const auto x_coord = left_end + static_cast(n) * dx; ftmp << x_coord << " " << Thyra::get_ele(x,n) << std::endl; } ftmp.close(); @@ -637,6 +624,36 @@ TEUCHOS_UNIT_TEST(BDF2, CDR) Teuchos::TimeMonitor::summarize(); } +// ************************************************************ +// ************************************************************ +TEUCHOS_UNIT_TEST(BDF2, CDR) +{ + // Create a communicator for Epetra objects + RCP comm; +#ifdef Tempus_ENABLE_MPI + comm = rcp(new Epetra_MpiComm(MPI_COMM_WORLD)); +#else + comm = rcp(new Epetra_SerialComm); +#endif + + CDR_Test>(comm, comm->NumProc(), out, success); +} + +// ************************************************************ +// ************************************************************ +TEUCHOS_UNIT_TEST(BDF2, CDR_Tpetra) +{ + // Get default Tpetra template types + using SC = Tpetra::Vector<>::scalar_type; + using LO = Tpetra::Vector<>::local_ordinal_type; + using GO = Tpetra::Vector<>::global_ordinal_type; + using Node = Tpetra::Vector<>::node_type; + + auto comm = Tpetra::getDefaultComm(); + + CDR_Test>(comm, comm->getSize(), out, success); +} + // ************************************************************ // ************************************************************