diff --git a/.github/workflows/codeql.yml b/.github/workflows/codeql.yml index 7176a77b0149..6e067a99b3bb 100644 --- a/.github/workflows/codeql.yml +++ b/.github/workflows/codeql.yml @@ -69,7 +69,7 @@ jobs: name: Get dependencies working-directory: ./packages/framework run: | - ./get_dependencies.sh --container + ./get_dependencies.sh - if: matrix.build-mode == 'manual' name: Generate CMake fragments diff --git a/packages/PyTrilinos2/CMakeLists.txt b/packages/PyTrilinos2/CMakeLists.txt index 3813a95d5cf0..3097975cfea5 100644 --- a/packages/PyTrilinos2/CMakeLists.txt +++ b/packages/PyTrilinos2/CMakeLists.txt @@ -1,6 +1,3 @@ -# -*- cmake -*- - - # Set the package name TRIBITS_PACKAGE(PyTrilinos2 DISABLE_STRONG_WARNINGS) @@ -27,11 +24,11 @@ PYTRILINOS2_CMAKE_ERROR ON ) TRIBITS_ADD_OPTION_AND_DEFINE(PyTrilinos2_BINDER_VERBOSE - PYTRILINOS2_B_VERBOSE + PYTRILINOS2_VERBOSE "Increase the verbosity of binder." OFF ) -SET(PyTrilinos2_BINDER_NUM_FILES "100" CACHE STRING "Maxinum number of generated files by binder.") +SET(PyTrilinos2_BINDER_NUM_FILES "150" CACHE STRING "Maxinum number of generated files by binder.") MESSAGE("-- Python3_EXECUTABLE:") IF(NOT DEFINED ${Python3_EXECUTABLE}) @@ -212,6 +209,10 @@ add_custom_command( add_custom_target(generate_include_name DEPENDS ${binder_include_name}) add_dependencies(generate_include_name generate_ETI_name) +ASSERT_DEFINED( + ${PACKAGE_NAME}_ENABLE_MueLu +) + set(BINDER_OPTIONS "") list(APPEND BINDER_OPTIONS --root-module PyTrilinos2) list(APPEND BINDER_OPTIONS --prefix ${CMAKE_CURRENT_BINARY_DIR}/binder) @@ -223,8 +224,13 @@ ELSE() ENDIF() list(APPEND BINDER_OPTIONS --bind Teuchos) list(APPEND BINDER_OPTIONS --bind Tpetra) -list(APPEND BINDER_OPTIONS --bind MueLu) -IF(PYTRILINOS2_B_VERBOSE) +list(APPEND BINDER_OPTIONS --bind Thyra) +list(APPEND BINDER_OPTIONS --bind ThyraTpetraAdapters) +list(APPEND BINDER_OPTIONS --bind Stratimikos) +IF(${PACKAGE_NAME}_ENABLE_MueLu) + list(APPEND BINDER_OPTIONS --bind MueLu) +ENDIF() +IF(PYTRILINOS2_VERBOSE) list(APPEND BINDER_OPTIONS -v) ENDIF() IF(PYTRILINOS2_SUPPRESS_ERRORS) @@ -346,7 +352,9 @@ ADD_SUBDIRECTORY( src ) file(COPY ${PyTrilinos2PyFiles} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/PyTrilinos2/.) #file(COPY ${PyTrilinos2PyFilesSo} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/PyTrilinos2/.) -SET(PyTrilinos2_PYTHONPATH "${CMAKE_CURRENT_BINARY_DIR}:$ENV{PYTHONPATH}") +SET(PyTrilinos2_PYTHONPATH "${CMAKE_CURRENT_BINARY_DIR}") + +TRIBITS_ADD_EXAMPLE_DIRECTORIES(examples) TRIBITS_ADD_TEST_DIRECTORIES(test) diff --git a/packages/PyTrilinos2/cmake/Dependencies.cmake b/packages/PyTrilinos2/cmake/Dependencies.cmake index 757550571c6c..9ef42d6ed3b6 100644 --- a/packages/PyTrilinos2/cmake/Dependencies.cmake +++ b/packages/PyTrilinos2/cmake/Dependencies.cmake @@ -10,9 +10,11 @@ SET(LIB_REQUIRED_DEP_PACKAGES Teuchos Tpetra - MueLu - ) -SET(LIB_OPTIONAL_DEP_PACKAGES) + Thyra + ThyraTpetraAdapters + Stratimikos +) +SET(LIB_OPTIONAL_DEP_PACKAGES MueLu) SET(TEST_REQUIRED_DEP_PACKAGES) SET(TEST_OPTIONAL_DEP_PACKAGES) SET(LIB_REQUIRED_DEP_TPLS) diff --git a/packages/PyTrilinos2/cmake/PyTrilinos2MakeTest.cmake b/packages/PyTrilinos2/cmake/PyTrilinos2MakeTest.cmake index 52691290a084..9b6c1e266e41 100644 --- a/packages/PyTrilinos2/cmake/PyTrilinos2MakeTest.cmake +++ b/packages/PyTrilinos2/cmake/PyTrilinos2MakeTest.cmake @@ -1,8 +1,6 @@ MACRO(PyTrilinos2_MAKE_MPI_TEST TEST_NAME) - FILE(COPY ${CMAKE_CURRENT_SOURCE_DIR}/${TEST_NAME}.py DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) - TRIBITS_ADD_TEST( ${Python3_EXECUTABLE} NOEXEPREFIX @@ -10,8 +8,12 @@ MACRO(PyTrilinos2_MAKE_MPI_TEST TEST_NAME) NAME ${TEST_NAME} ARGS "${TEST_NAME}.py" PASS_REGULAR_EXPRESSION "OK" - ENVIRONMENT "PYTHONPATH=${PyTrilinos2_PYTHONPATH}" + ENVIRONMENT "PYTHONPATH=${PyTrilinos2_PYTHONPATH}:$ENV{PYTHONPATH}" ${ARGN} ) + TRIBITS_COPY_FILES_TO_BINARY_DIR(${TEST_NAME}_cp + SOURCE_FILES ${TEST_NAME}.py + CATEGORIES BASIC) + ENDMACRO(PyTrilinos2_MAKE_MPI_TEST TEST_NAME) diff --git a/packages/PyTrilinos2/cmake/PyTrilinos2_config.hpp.in b/packages/PyTrilinos2/cmake/PyTrilinos2_config.hpp.in new file mode 100644 index 000000000000..1f08c770a414 --- /dev/null +++ b/packages/PyTrilinos2/cmake/PyTrilinos2_config.hpp.in @@ -0,0 +1,6 @@ +#ifndef PYTILINOS2_CONFIG_HPP +#define PYTILINOS2_CONFIG_HPP + +#cmakedefine HAVE_PYTRILINOS2_MUELU + +#endif diff --git a/packages/PyTrilinos2/examples/CMakeLists.txt b/packages/PyTrilinos2/examples/CMakeLists.txt new file mode 100644 index 000000000000..b6ec0418afb4 --- /dev/null +++ b/packages/PyTrilinos2/examples/CMakeLists.txt @@ -0,0 +1,3 @@ +INCLUDE(PyTrilinos2MakeTest) + +PyTrilinos2_MAKE_MPI_TEST(exampleCG) diff --git a/packages/PyTrilinos2/examples/CG.py b/packages/PyTrilinos2/examples/exampleCG.py similarity index 97% rename from packages/PyTrilinos2/examples/CG.py rename to packages/PyTrilinos2/examples/exampleCG.py index e968eff905ef..845ae8d9971a 100644 --- a/packages/PyTrilinos2/examples/CG.py +++ b/packages/PyTrilinos2/examples/exampleCG.py @@ -149,6 +149,14 @@ def main(): plt.plot(x0_view) plt.savefig('x0_view.png', dpi=800, bbox_inches='tight',pad_inches = 0) + success = True + if comm.getRank() == 0: + if success: + print("OK") + else: + print("FAIL") + + if __name__ == "__main__": comm = MPI.COMM_WORLD rank = comm.Get_rank() diff --git a/packages/PyTrilinos2/python/__init__.py b/packages/PyTrilinos2/python/__init__.py index 696783bd3b7c..a9418bf17c8b 100644 --- a/packages/PyTrilinos2/python/__init__.py +++ b/packages/PyTrilinos2/python/__init__.py @@ -10,3 +10,5 @@ __version__ = '0.1.0' def version(): return 'PyTrilinos2 version: ' + __version__ + +from . PyTrilinos2 import * diff --git a/packages/PyTrilinos2/scripts/PyTrilinos2_config.cfg b/packages/PyTrilinos2/scripts/PyTrilinos2_config.cfg index 148092c485a9..185f0bb8bbb2 100644 --- a/packages/PyTrilinos2/scripts/PyTrilinos2_config.cfg +++ b/packages/PyTrilinos2/scripts/PyTrilinos2_config.cfg @@ -9,23 +9,77 @@ +include +include + +# Standard library +-class std::basic_ios +-class std::vector +-class std::map +-class std::integral_constant +-class std::iterator +-class std::reverse_iterator +-class std::set + +# OpenMPI +-class ompi_status_public_t +-class ompi_request_t +-class ompi_errhandler_t +-class ompi_communicator_t + +############################################################################### +# Teuchos +include -+custom_shared Teuchos::RCP + +include_for_namespace Teuchos +add_on_binder_for_namespace Teuchos def_Teuchos_functions -+add_on_binder_for_namespace Tpetra def_initialize_Kokkos + ++custom_shared Teuchos::RCP + +-namespace Teuchos::Details + +-function Teuchos::Ptr::operator-> +-function Teuchos::Ptr::get +-function Teuchos::Ptr::getRawPtr +-function Teuchos::Array::data +-function Teuchos::Array::getRawPtr +-function Teuchos::ArrayRCP::getRawPtr +-function Teuchos::ArrayView::access_private_ptr +-function Teuchos::ArrayView::begin +-function Teuchos::ArrayView::end +-function Teuchos::ArrayView::data +-function Teuchos::ArrayView::getRawPtr +-class Teuchos::ArrayRCP +-class Teuchos::ArrayRCP +-class Teuchos::Array +-class Teuchos::Ptr +-class Teuchos::Ptr +-class Teuchos::Ptr +-class Teuchos::Ptr +-class Teuchos::Ptr +-class Teuchos::Ptr +-class Teuchos::Ptr +-class Teuchos::Ptr +-class Teuchos::Ptr + ++add_on_binder Teuchos::ParameterList def_ParameterList_member_functions -function Teuchos::ParameterList::sublist -function Teuchos::ParameterList::set -function Teuchos::ParameterList::get -+add_on_binder Teuchos::ParameterList def_ParameterList_member_functions -+include_for_namespace Tpetra -+add_on_binder Tpetra::CrsGraph define_CrsGraph_member_functions -+include -+add_on_binder Tpetra::CrsMatrix define_CrsMatrix_member_functions -+add_on_binder Tpetra::Vector define_Vector_member_functions -+add_on_binder Tpetra::MultiVector define_MultiVector_member_functions --class Kokkos::View --class Kokkos::DualView +-class Teuchos::ParameterListAcceptorDefaultBase + +-class Teuchos::PromotionTraits +-function Teuchos::Serializer::createObj +-function Teuchos::TimeMonitor::computeGlobalTimerStatistics +-class Teuchos::CommandLineProcessor::enum_opt_data_t +-class Teuchos::CommandLineProcessor::TimeMonitorSurrogate +-class Teuchos::RawWorkspace +-function Teuchos::getRawMpiComm +-function Teuchos::rcp_dynamic_cast +-function Teuchos::ptr_dynamic_cast +-class Teuchos::BLAS +-class Teuchos::Describable + +############################################################################### +# Kokkos -namespace Kokkos::Impl -namespace KokkosBlas -include @@ -91,88 +145,71 @@ -include -include -include --class Teuchos::ArrayView --class Teuchos::ArrayView --class Teuchos::Serializer --class Teuchos::PromotionTraits --class Teuchos::CommandLineProcessor --class Teuchos::Ptr >> --class Teuchos::ArrayView >> --class Teuchos::TimeMonitor +-include +-include +-include +-class Kokkos::Device +-class Kokkos::DualView +-class Kokkos::DynRankView +-class Kokkos::View + +############################################################################### +# Tpetra ++add_on_binder_for_namespace Tpetra def_initialize_Kokkos ++include_for_namespace Tpetra ++include ++add_on_binder Tpetra::CrsGraph define_CrsGraph_member_functions ++add_on_binder Tpetra::Vector define_Vector_member_functions ++add_on_binder Tpetra::MultiVector define_MultiVector_member_functions ++add_on_binder Tpetra::CrsMatrix define_CrsMatrix_member_functions -namespace Tpetra::Details -namespace Tpetra::Import_Util --namespace Teuchos::Details -namespace Tpetra::KokkosRefactor --function Tpetra::Details::isInterComm --function Tpetra::SrcDistObject::operator= --function Teuchos::TimeMonitor::computeGlobalTimerStatistics --function Teuchos::mpiErrorCodeToString --class Teuchos::CommandLineProcessor::enum_opt_data_t --class Teuchos::CommandLineProcessor::TimeMonitorSurrogate --class Teuchos::RawWorkspace --class std::ostream --class std::basic_ios --class std::vector --class std::map --class std::integral_constant --class std::integral_constant --class std::iterator --class std::reverse_iterator --class Teuchos::ArrayView --class Teuchos::ArrayRCP --class Teuchos::Array --class Teuchos::Describable --class Teuchos::BLAS --function Teuchos::fancyOStream --class Teuchos::ParameterListAcceptorDefaultBase --class Teuchos::Dependency --class Teuchos::DependencySheet --class Teuchos::MpiCommRequestBase --function Teuchos::getRawMpiComm --class Teuchos::OpaqueWrapper --class Teuchos::OpaqueWrapper --function Teuchos::Details::setMpiReductionOp --function Teuchos::Details::getMpiOpForEReductionType --function Tpetra::Details::extractMpiCommFromTeuchos --class Teuchos::OpaqueWrapper --function Teuchos::Details::setMpiReductionOp --function Teuchos::Details::getMpiOpForEReductionType --function Tpetra::Details::PackCrsGraphImpl::packRow +-class Tpetra::BlockCrsMatrix -class Tpetra::Details::DistributorPlan +-class Tpetra::Directory +-class Tpetra::Distribution +-class Tpetra::Distribution1D +-class Tpetra::Distribution2D +-class Tpetra::DistributionLowerTriangularBlock +-class Tpetra::DistributionMM -class Tpetra::DistributionType --namespace Xpetra --class Thyra::VectorSpaceBase --class Thyra::ProductVectorSpace --class Tpetra::Details::DistributorPlan -class Tpetra::Distributor --class Tpetra::DistributionType +-class Tpetra::LowerTriangularBlockOperator -class Tpetra::distributorSendTypes +-class Tpetra::ImportExportData +-function Tpetra::Details::extractMpiCommFromTeuchos +-function Tpetra::Details::isInterComm +-function Tpetra::Details::PackCrsGraphImpl::packRow +-function Tpetra::SrcDistObject::operator= + +############################################################################### +# Xpetra +-namespace Xpetra + +############################################################################### +# MueLu +-class MueLu::BaseClass +-class MueLu::Describable -class MueLu::FactoryAcceptor +-class MueLu::FactoryBase -class MueLu::FactoryFactory --class MueLu::FactoryManagerBase -class MueLu::FactoryManager --class MueLu::FactoryBase +-class MueLu::FactoryManagerBase -class MueLu::Hierarchy -class MueLu::HierarchyManager +-class MueLu::Level -class MueLu::TimeMonitor --class MueLu::Describable --class Kokkos::Device --class Tpetra::DistributionType --class Tpetra::Directory --class Tpetra::Distribution --class Tpetra::Distribution1D --class Tpetra::Distribution2D --class Tpetra::DistributionMM --class Tpetra::DistributionLowerTriangularBlock --class MueLu::BaseClass --function Teuchos::rcp_dynamic_cast --class Teuchos::VerboseObjectBase --class MueLu::VerboseObject --class Tpetra::BlockCrsMatrix --class Kokkos::DynRankView -class MueLu::VariableContainer --class Tpetra::LowerTriangularBlockOperator --class MueLu::Level --class MueLu::Hierarchy --class std::set --class Teuchos::Ptr +-class MueLu::VerboseObject ++include + +############################################################################### +# Thyra ++include +-class Thyra::ProductVectorSpace +-class Thyra::RowStatLinearOpBase +-class Thyra::ScaledLinearOpBase +-class Thyra::LinearOpSourceBase ++include_for_namespace Thyra ++add_on_binder Thyra::LinearOpWithSolveBase define_solve \ No newline at end of file diff --git a/packages/PyTrilinos2/src/CMakeLists.txt b/packages/PyTrilinos2/src/CMakeLists.txt index 712db1dc6a96..f9931f40104d 100644 --- a/packages/PyTrilinos2/src/CMakeLists.txt +++ b/packages/PyTrilinos2/src/CMakeLists.txt @@ -1,6 +1,7 @@ +TRIBITS_CONFIGURE_FILE(${PACKAGE_NAME}_config.hpp) FILE(GLOB PYTRILINOS2_SRC "${CMAKE_CURRENT_SOURCE_DIR}/*.cpp") -MESSAGE("PYTRILINOS2_SRC = ${PYTRILINOS2_SRC}") +# MESSAGE("PYTRILINOS2_SRC = ${PYTRILINOS2_SRC}") FILE(COPY ${PYTRILINOS2_SRC} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) list(APPEND PYTRILINOS2_SRC ${CMAKE_CURRENT_BINARY_DIR}/../binder/PyTrilinos2.cpp) @@ -23,7 +24,7 @@ IF(NOT PYTRILINOS2_USE_ONE_FILE) ) ENDIF() -MESSAGE("PYTRILINOS2_SRC with binder = ${PYTRILINOS2_SRC}") +# MESSAGE("PYTRILINOS2_SRC with binder = ${PYTRILINOS2_SRC}") pybind11_add_module(PyTrilinos2 ${PYTRILINOS2_SRC}) @@ -31,7 +32,7 @@ add_dependencies(PyTrilinos2 binder_call generate_ETI_name generate_include_name target_include_directories(PyTrilinos2 PUBLIC ${Mpi4Py_INCLUDE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}) target_compile_features(PyTrilinos2 PUBLIC cxx_std_14) -foreach(depPkg IN LISTS PyTrilinos2_LIB_ENABLED_DEPENDENCIES) +foreach(depPkg IN LISTS PyTrilinos2_LIB_ENABLED_DEPENDENCIES) target_link_libraries(PyTrilinos2 PUBLIC ${depPkg}::all_libs) endforeach() target_link_libraries(PyTrilinos2 PUBLIC ${Trilinos_EXTRA_LINK_FLAGS}) @@ -43,4 +44,4 @@ INSTALL(TARGETS PyTrilinos2 add_custom_command(TARGET PyTrilinos2 POST_BUILD COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_BINARY_DIR}/PyTrilinos2.so ${CMAKE_CURRENT_BINARY_DIR}/../PyTrilinos2/. COMMENT "Copy ${PROJECT_BINARY_DIR}/src/PyTrilinos2.so" -) \ No newline at end of file +) diff --git a/packages/PyTrilinos2/src/PyTrilinos2_Binder_Input.hpp b/packages/PyTrilinos2/src/PyTrilinos2_Binder_Input.hpp index 0b313765deee..a82036562be6 100644 --- a/packages/PyTrilinos2/src/PyTrilinos2_Binder_Input.hpp +++ b/packages/PyTrilinos2/src/PyTrilinos2_Binder_Input.hpp @@ -7,6 +7,12 @@ // ***************************************************************************** // @HEADER +#include #include #include +#include +#include +#ifdef HAVE_PYTRILINOS2_MUELU #include +#endif +#include diff --git a/packages/PyTrilinos2/src/PyTrilinos2_Stratimikos_ETI.hpp b/packages/PyTrilinos2/src/PyTrilinos2_Stratimikos_ETI.hpp new file mode 100644 index 000000000000..ee779472d873 --- /dev/null +++ b/packages/PyTrilinos2/src/PyTrilinos2_Stratimikos_ETI.hpp @@ -0,0 +1,42 @@ +// @HEADER +// ***************************************************************************** +// PyTrilinos2: Automatic Python Interfaces to Trilinos Packages +// +// Copyright 2022 NTESS and the PyTrilinos2 contributors. +// SPDX-License-Identifier: BSD-3-Clause +// ***************************************************************************** +// @HEADER + +#ifndef PYTRILINOS2_STRATIMIKOS_ETI +#define PYTRILINOS2_STRATIMIKOS_ETI + +#include +#include +#ifdef HAVE_PYTRILINOS2_MUELU +#include +#endif + + +#define BINDER_STRATIMIKOS_LINEARSOLVERBUILDER_INSTANT(SC) \ + template class Stratimikos::LinearSolverBuilder; + +#define BINDER_STRATIMIKOS_MUELU_INSTANT(SCALAR,LO,GO,NO) \ + template void enableMueLu(Stratimikos::LinearSolverBuilder& builder, const std::string& stratName); \ + template void enableMueLuRefMaxwell(Stratimikos::LinearSolverBuilder& builder, const std::string& stratName); \ + template void enableMueLuMaxwell1(Stratimikos::LinearSolverBuilder& builder, const std::string& stratName); + + +namespace Stratimikos { + + template + void initiate(T) {}; + + BINDER_STRATIMIKOS_LINEARSOLVERBUILDER_INSTANT(Tpetra::Details::DefaultTypes::scalar_type) + +#ifdef HAVE_PYTRILINOS2_MUELU + BINDER_STRATIMIKOS_MUELU_INSTANT(Tpetra::Details::DefaultTypes::scalar_type, Tpetra::Details::DefaultTypes::local_ordinal_type, Tpetra::Details::DefaultTypes::global_ordinal_type, Tpetra::KokkosClassic::DefaultNode::DefaultNodeType) +#endif + +} + +#endif // PYTRILINOS2_STRATIMIKOS_ETI diff --git a/packages/PyTrilinos2/src/PyTrilinos2_Thyra_Custom.hpp b/packages/PyTrilinos2/src/PyTrilinos2_Thyra_Custom.hpp new file mode 100644 index 000000000000..e7e1603c1eba --- /dev/null +++ b/packages/PyTrilinos2/src/PyTrilinos2_Thyra_Custom.hpp @@ -0,0 +1,17 @@ +#ifndef PYTRILINOS2_THYRA_CUSTOM_HPP +#define PYTRILINOS2_THYRA_CUSTOM_HPP + +#include + +template +void define_solve(T cl) { + using SCALAR = typename T::type::scalar_type; + + cl.def("solve",[](Teuchos::RCP > &m, + const Thyra::EOpTransp A_trans, + const Thyra::MultiVectorBase &B, + const Teuchos::RCP > &X) + { return m->solve(A_trans, B, X.ptr()); }); +} + +#endif diff --git a/packages/PyTrilinos2/src/PyTrilinos2_Thyra_ETI.hpp b/packages/PyTrilinos2/src/PyTrilinos2_Thyra_ETI.hpp new file mode 100644 index 000000000000..4de39cb88d45 --- /dev/null +++ b/packages/PyTrilinos2/src/PyTrilinos2_Thyra_ETI.hpp @@ -0,0 +1,87 @@ +// @HEADER +// ***************************************************************************** +// PyTrilinos2: Automatic Python Interfaces to Trilinos Packages +// +// Copyright 2022 NTESS and the PyTrilinos2 contributors. +// SPDX-License-Identifier: BSD-3-Clause +// ***************************************************************************** +// @HEADER + +#ifndef PYTRILINOS2_THYRA_ETI +#define PYTRILINOS2_THYRA_ETI + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + + +#define BINDER_THYRA_WRAPPERS_INSTANT(SC,LO,GO,NO) \ + template Teuchos::RCP > createVector(const Teuchos::RCP > &tpetraVector, \ + const Teuchos::RCP > space); \ + template Teuchos::RCP > createMultiVector(const Teuchos::RCP > &tpetraMultiVector, \ + const Teuchos::RCP > rangeSpace, \ + const Teuchos::RCP > domainSpace); + +#define BINDER_THYRA_EUCLIDEANSCALARPROD_INSTANT(SC,LO,GO,NO) \ + template class Thyra::TpetraEuclideanScalarProd; + +#define BINDER_THYRA_VECTORSPACE_INSTANT(SC,LO,GO,NO) \ + template class Thyra::TpetraVectorSpace; + +#define BINDER_THYRA_MULTIVECTOR_INSTANT(SC,LO,GO,NO) \ + template class Thyra::TpetraMultiVector; + +#define BINDER_THYRA_VECTOR_INSTANT(SC, LO, GO, NO) \ + template class Thyra::TpetraVector; + +#define BINDER_THYRA_LINEAROP_INSTANT(SC,LO,GO,NO) \ + template class Thyra::TpetraLinearOp; \ + template Teuchos::RCP > tpetraLinearOp(const Teuchos::RCP > &rangeSpace, \ + const Teuchos::RCP > &domainSpace, \ + const Teuchos::RCP > &tpetraOperator); \ + template Teuchos::RCP > constTpetraLinearOp(const Teuchos::RCP > &rangeSpace, \ + const Teuchos::RCP > &domainSpace, \ + const Teuchos::RCP > &tpetraOperator); + +#define BINDER_THYRA_LINEARSOLVER_INSTANT(SC) \ + template class Thyra::LinearOpWithSolveFactoryBase; \ + template Teuchos::RCP > Thyra::createLinearSolveStrategy(const Thyra::LinearSolverBuilderBase &linearSolverBuilder, const std::string &linearSolveStrategyName); \ + template Teuchos::RCP > linearOpWithSolve(const Thyra::LinearOpWithSolveFactoryBase &lowsFactory, \ + const Teuchos::RCP > &fwdOp, \ + const Thyra::ESupportSolveUse supportSolveUse); \ + template struct Thyra::SolveCriteria; \ + template struct Thyra::SolveStatus; + +#define BINDER_TPETRA_OPERATOR_INSTANT(SC,LO,GO,NO) \ + template class Tpetra::Operator; + +namespace Tpetra { + BINDER_TPETRA_OPERATOR_INSTANT(Tpetra::Details::DefaultTypes::scalar_type, Tpetra::Details::DefaultTypes::local_ordinal_type, Tpetra::Details::DefaultTypes::global_ordinal_type, Tpetra::Details::DefaultTypes::node_type) +} + + +namespace Thyra { + + template + void initiate(T) {}; + + BINDER_THYRA_WRAPPERS_INSTANT(Tpetra::Details::DefaultTypes::scalar_type, Tpetra::Details::DefaultTypes::local_ordinal_type, Tpetra::Details::DefaultTypes::global_ordinal_type, Tpetra::Details::DefaultTypes::node_type) + BINDER_THYRA_EUCLIDEANSCALARPROD_INSTANT(Tpetra::Details::DefaultTypes::scalar_type, Tpetra::Details::DefaultTypes::local_ordinal_type, Tpetra::Details::DefaultTypes::global_ordinal_type, Tpetra::Details::DefaultTypes::node_type) + BINDER_THYRA_VECTORSPACE_INSTANT(Tpetra::Details::DefaultTypes::scalar_type, Tpetra::Details::DefaultTypes::local_ordinal_type, Tpetra::Details::DefaultTypes::global_ordinal_type, Tpetra::Details::DefaultTypes::node_type) + BINDER_THYRA_MULTIVECTOR_INSTANT(Tpetra::Details::DefaultTypes::scalar_type, Tpetra::Details::DefaultTypes::local_ordinal_type, Tpetra::Details::DefaultTypes::global_ordinal_type, Tpetra::Details::DefaultTypes::node_type) + BINDER_THYRA_VECTOR_INSTANT(Tpetra::Details::DefaultTypes::scalar_type, Tpetra::Details::DefaultTypes::local_ordinal_type, Tpetra::Details::DefaultTypes::global_ordinal_type, Tpetra::Details::DefaultTypes::node_type) + BINDER_THYRA_LINEAROP_INSTANT(Tpetra::Details::DefaultTypes::scalar_type, Tpetra::Details::DefaultTypes::local_ordinal_type, Tpetra::Details::DefaultTypes::global_ordinal_type, Tpetra::Details::DefaultTypes::node_type) + BINDER_THYRA_LINEARSOLVER_INSTANT(Tpetra::Details::DefaultTypes::scalar_type) +} + +#endif // PYTRILINOS2_THYRA_ETI diff --git a/packages/PyTrilinos2/src/PyTrilinos2_Tpetra_Custom.hpp b/packages/PyTrilinos2/src/PyTrilinos2_Tpetra_Custom.hpp index 0cc6458aac2d..a24fb33efd96 100644 --- a/packages/PyTrilinos2/src/PyTrilinos2_Tpetra_Custom.hpp +++ b/packages/PyTrilinos2/src/PyTrilinos2_Tpetra_Custom.hpp @@ -316,18 +316,15 @@ void define_MultiVector_member_functions(T cl) { template void def_initialize_Kokkos(T m) { m.def("initialize_Kokkos",[](int num_threads, - int num_devices, int device_id){ if(!Kokkos::is_initialized()) { Kokkos::InitializationSettings args; args.set_num_threads(num_threads); - args.set_num_devices(num_devices); args.set_device_id(device_id); Kokkos::initialize(args); } }, py::arg("num_threads") = -1, - py::arg("num_devices") = -1, py::arg("device_id") = -1 ); m.def("finalize_Kokkos",[](){ diff --git a/packages/PyTrilinos2/src/PyTrilinos2_stream.hpp b/packages/PyTrilinos2/src/PyTrilinos2_stream.hpp new file mode 100644 index 000000000000..7ff146f91798 --- /dev/null +++ b/packages/PyTrilinos2/src/PyTrilinos2_stream.hpp @@ -0,0 +1,26 @@ +#ifndef TEUCHOS_STREAM_HPP +#define TEUCHOS_STREAM_HPP + +#include +#include +#include "Teuchos_RCP.hpp" + +namespace Teuchos { + + std::ofstream openOfstream(std::string filename){ + std::ofstream outdata; + outdata.open(filename); + return outdata; + } + + void closeOfstream(std::ofstream &outdata) { + outdata.close(); + } + + Teuchos::RCP getCout() { + Teuchos::RCP out = Teuchos::rcp(&std::cout, false); + return out; + } +} + +#endif diff --git a/packages/PyTrilinos2/test/CMakeLists.txt b/packages/PyTrilinos2/test/CMakeLists.txt index 385477db14c1..b65b8fd9bd18 100644 --- a/packages/PyTrilinos2/test/CMakeLists.txt +++ b/packages/PyTrilinos2/test/CMakeLists.txt @@ -5,3 +5,21 @@ INCLUDE(PyTrilinos2MakeTest) PyTrilinos2_MAKE_MPI_TEST(CG) PyTrilinos2_MAKE_MPI_TEST(parameterList) + +TRIBITS_COPY_FILES_TO_BINARY_DIR(Stratimiko_cp + SOURCE_FILES Stratimikos.py) + +TRIBITS_ADD_TEST( + ${Python3_EXECUTABLE} + NOEXEPREFIX + NOEXESUFFIX + NAME Stratimikos + POSTFIX_AND_ARGS_0 "LU" Stratimikos.py --problemSize=100 --solver=LU + POSTFIX_AND_ARGS_1 "CG" Stratimikos.py --problemSize=100 --solver=CG --prec=None + POSTFIX_AND_ARGS_2 "CG_Jacobi" Stratimikos.py --problemSize=100 --solver=CG --prec=Jacobi + POSTFIX_AND_ARGS_3 "BiCGStab_Chebyshev" Stratimikos.py --problemSize=100 --solver=BiCGStab --prec=Chebyshev + POSTFIX_AND_ARGS_4 "GMRES_ILU" Stratimikos.py --problemSize=1000 --solver=GMRES --prec=ILU + POSTFIX_AND_ARGS_5 "CG_multigrid" Stratimikos.py --problemSize=10000 --solver=CG --prec=multigrid + PASS_REGULAR_EXPRESSION "OK" + ENVIRONMENT "PYTHONPATH=${PyTrilinos2_PYTHONPATH}:$ENV{PYTHONPATH}" +) diff --git a/packages/PyTrilinos2/test/Stratimikos.py b/packages/PyTrilinos2/test/Stratimikos.py new file mode 100755 index 000000000000..7486846ceacf --- /dev/null +++ b/packages/PyTrilinos2/test/Stratimikos.py @@ -0,0 +1,253 @@ +#!/usr/env python +# @HEADER +# ***************************************************************************** +# PyTrilinos2: Automatic Python Interfaces to Trilinos Packages +# +# Copyright 2022 NTESS and the PyTrilinos2 contributors. +# SPDX-License-Identifier: BSD-3-Clause +# ***************************************************************************** +# @HEADER + +from mpi4py import MPI +from argparse import ArgumentParser +try: + # MPI, Timers, ParameterList + from PyTrilinos2 import Teuchos + # Linear algebra + from PyTrilinos2 import Tpetra + # Linear algebra interfaces used by Stratimikos + from PyTrilinos2 import Thyra + # Unified solver & preconditioner interface + from PyTrilinos2 import Stratimikos + from PyTrilinos2.getTpetraTypeName import getTypeName, getDefaultNodeType +except ImportError: + print("\nFailed to import PyTrilinos2. Consider setting the Python load path in your environment with\n export PYTHONPATH=${TRILINOS_BUILD_DIR}/packages/PyTrilinos2:${PYTHONPATH}\nwhere TRILINOS_BUILD_DIR is the build directory of your Trilinos build.\n") + raise + +try: + import matplotlib as mpl + mpl.use('Agg') + mpl.rcParams.update(mpl.rcParamsDefault) + import matplotlib.pyplot as plt + display = True +except: + display = False + + +def assemble1DLaplacian(globalNumRows, comm): + Tpetra_Map = getTypeName('Map') + Tpetra_CrsGraph = getTypeName('CrsGraph') + Tpetra_CrsMatrix = getTypeName('CrsMatrix') + + rowmap = Tpetra_Map(globalNumRows, 0, comm) + + graph = Tpetra_CrsGraph(rowmap, 3) + for lid in range(rowmap.getMinLocalIndex(), rowmap.getMaxLocalIndex()+1): + gid = rowmap.getGlobalElement(lid) + indices = [gid] + if gid > 0: + indices.append(gid-1) + if gid < rowmap.getMaxAllGlobalIndex(): + indices.append(gid+1) + graph.insertGlobalIndices(gid, indices) + graph.fillComplete() + + A = Tpetra_CrsMatrix(graph) + for lid in range(rowmap.getMinLocalIndex(), rowmap.getMaxLocalIndex()+1): + gid = rowmap.getGlobalElement(lid) + indices = [gid] + vals = [2.] + if gid > 0: + indices.append(gid-1) + vals.append(-1.) + if gid < rowmap.getMaxAllGlobalIndex(): + indices.append(gid+1) + vals.append(-1.) + A.replaceGlobalValues(gid, indices, vals) + A.fillComplete() + return A + + +def main(): + comm = Teuchos.getTeuchosComm(MPI.COMM_WORLD) + + parser = ArgumentParser() + parser.add_argument("--problemSize", default=1000, help="Global problem size", type=int) + parser.add_argument("--solver", default="GMRES", choices=["LU", "CG", "BiCGStab", "GMRES"], help="Linear solver") + parser.add_argument("--prec", default="None", choices=["None", "Jacobi", "Chebyshev", "ILU", "multigrid"], help="Preconditioner") + parser.add_argument("--maxIts", default=100, help="Maximum number of iterations", type=int) + args = parser.parse_args() + + # Set up timers + timer = Teuchos.StackedTimer("Main") + Teuchos.TimeMonitor.setStackedTimer(timer) + + # Set up output streams + cout = Teuchos.getCout() + out = Teuchos.fancyOStream(cout) + + # Set up Tpetra matrices and vectors + Tpetra_Map = getTypeName('Map') + Tpetra_Vector = getTypeName('Vector') + Tpetra_Export = getTypeName('Export') + + tpetra_A = assemble1DLaplacian(args.problemSize, comm) + tpetra_map = tpetra_A.getRowMap() + + n0 = (args.problemSize if comm.getRank() == 0 else 0) + tpetra_map0 = Tpetra_Map(args.problemSize, n0, 0, comm) + + tpetra_x = Tpetra_Vector(tpetra_map, True) + tpetra_b = Tpetra_Vector(tpetra_map, False) + tpetra_residual = Tpetra_Vector(tpetra_map, False) + + tpetra_b.putScalar(1.) + + tpetra_A.describe(out) + + # Wrap Tpetra objects as Thyra objects + thyra_map = Thyra.tpetraVectorSpace(tpetra_map) + thyra_x = Thyra.tpetraVector(thyra_map, tpetra_x) + thyra_b = Thyra.tpetraVector(thyra_map, tpetra_b) + thyra_residual = Thyra.tpetraVector(thyra_map, tpetra_residual) + thyra_A = Thyra.tpetraLinearOp(thyra_map, thyra_map, tpetra_A) + + # Set up linear solver + linearSolverBuilder = Stratimikos.LinearSolverBuilder_double_t() + + # Hook up preconditioners that are not enabled by default + if hasattr(Stratimikos, "enableMueLu"): + Stratimikos.enableMueLu(linearSolverBuilder, "MueLu") + else: + assert args.prec != "multigrid", "\"multigrid\" preconditioner requires the MueLu package." + if hasattr(Stratimikos, "enableMueLuRefMaxwell"): + Stratimikos.enableMueLuRefMaxwell(linearSolverBuilder, "MueLuRefMaxwell") + if hasattr(Stratimikos, "enableMueLuMaxwell1"): + Stratimikos.enableMueLuMaxwell1(linearSolverBuilder, "MueLuMaxwell1") + + # Print default parameters + validParams = linearSolverBuilder.getValidParameters() + if comm.getRank() == 0: + print(validParams) + + params = Teuchos.ParameterList() + + # Set parameters for solver + if args.solver == "LU": + params["Linear Solver Type"] = "Amesos2" + elif args.solver in ("CG", "GMRES", "BiCGStab"): + if args.solver == "CG": + BelosSolver = "Pseudo Block CG" + elif args.solver == "GMRES": + BelosSolver = "Block GMRES" + elif args.solver == "BiCGStab": + BelosSolver = "BiCGStab" + params["Linear Solver Type"] = "Belos" + params["Linear Solver Types"] = Teuchos.ParameterList() + params["Linear Solver Types"]["Belos"] = Teuchos.ParameterList() + params["Linear Solver Types"]["Belos"]["Solver Type"] = BelosSolver + params["Linear Solver Types"]["Belos"]["Solver Types"] = Teuchos.ParameterList() + params["Linear Solver Types"]["Belos"]["Solver Types"][BelosSolver] = Teuchos.ParameterList() + params["Linear Solver Types"]["Belos"]["Solver Types"][BelosSolver]["Convergence Tolerance"] = 1e-8 + params["Linear Solver Types"]["Belos"]["Solver Types"][BelosSolver]["Maximum Iterations"] = args.maxIts + params["Linear Solver Types"]["Belos"]["Solver Types"][BelosSolver]["Verbosity"] = 41 + params["Linear Solver Types"]["Belos"]["Solver Types"][BelosSolver]["Output Frequency"] = 1 + params["Linear Solver Types"]["Belos"]["Solver Types"][BelosSolver]["Output Style"] = 1 + + params["Linear Solver Types"]["Belos"]["VerboseObject"] = Teuchos.ParameterList() + params["Linear Solver Types"]["Belos"]["VerboseObject"]["Verbosity Level"] = "low" + else: + raise NotImplementedError("Unknown solver: {}".format(args.solver)) + + # Set parameters for preconditioner + if args.prec == "None": + params["Preconditioner Type"] = "None" + elif args.prec in ("Jacobi", "ILU", "Chebyshev"): + params["Preconditioner Type"] = "Ifpack2" + params["Preconditioner Types"] = Teuchos.ParameterList() + params["Preconditioner Types"]["Ifpack2"] = Teuchos.ParameterList() + if args.prec == "Jacobi": + params["Preconditioner Types"]["Ifpack2"]["Prec Type"] = "relaxation" + params["Preconditioner Types"]["Ifpack2"]["Ifpack2 Settings"] = Teuchos.ParameterList() + params["Preconditioner Types"]["Ifpack2"]["Ifpack2 Settings"]["relaxation: type"] = "Jacobi" + params["Preconditioner Types"]["Ifpack2"]["Ifpack2 Settings"]["relaxation: sweeps"] = 1 + elif args.prec == "Chebyshev": + params["Preconditioner Types"]["Ifpack2"]["Prec Type"] = "chebyshev" + params["Preconditioner Types"]["Ifpack2"]["Ifpack2 Settings"] = Teuchos.ParameterList() + params["Preconditioner Types"]["Ifpack2"]["Ifpack2 Settings"]["chebyshev: degree"] = 2 + elif args.prec == "ILU": + params["Preconditioner Types"]["Ifpack2"]["Prec Type"] = "ILUT" + elif args.prec == "multigrid": + params["Preconditioner Type"] = "MueLu" + else: + raise NotImplementedError("Unknown preconditioner: {}".format(args.prec)) + + linearSolverBuilder.setParameterList(params) + + solverName = linearSolverBuilder.getLinearSolveStrategyName() + precName = linearSolverBuilder.getPreconditionerStrategyName() + solverFactory = Thyra.createLinearSolveStrategy(linearSolverBuilder) + + timer.start("Setup solver") + thyra_invA = Thyra.linearOpWithSolve(solverFactory, thyra_A) + timer.stop("Setup solver") + assert thyra_invA.solveSupports(Thyra.NOTRANS) + + # Solve + timer.start("Solve") + status = thyra_invA.solve(Thyra.NOTRANS, thyra_b, thyra_x) + timer.stop("Solve") + + # Compute residual + tpetra_A.apply(tpetra_x, tpetra_residual) + tpetra_residual.update(1, tpetra_b, -1) + resNorm = tpetra_residual.norm2() + + if comm.getRank() == 0: + print("Solver choice: ", args.solver) + print("Stratimikos solver name: ", solverName) + print("Preconditioner choice: ", args.prec) + print("Stratimikos preconditioner name: ", precName) + if solverName == "Belos": + its = status.extraParameters["Iteration Count"] + if comm.getRank() == 0: + print('Norm of residual after {} iterations = {} '.format(its, resNorm)) + elif comm.getRank() == 0: + print('Norm of residual = {} '.format(resNorm)) + + # Gather solution vector on rank 0 + tpetra_x0 = Tpetra_Vector(tpetra_map0, True) + export = Tpetra_Export(tpetra_map0, tpetra_map) + tpetra_x0.doImport(source=tpetra_x, exporter=export, CM=Tpetra.CombineMode.REPLACE) + + # Print timings + timer.stop("Main") + options = Teuchos.StackedTimer.OutputOptions() + options.output_fraction = options.output_histogram = options.output_minmax = True; + timer.report(out, comm, options) + + if comm.getRank() == 0 and display: + x0_view = tpetra_x0.getLocalViewHost() + plt.figure() + plt.plot(x0_view) + plt.savefig('x0_view.png', dpi=800, bbox_inches='tight', pad_inches=0) + + comm.barrier() + success = ((status.solveStatus == Thyra.ESolveStatus.SOLVE_STATUS_CONVERGED) and (resNorm < 1e-6)) + if comm.getRank() == 0: + if success: + print("OK") + else: + print("FAIL") + + +if __name__ == "__main__": + # initialize kokkos + if getDefaultNodeType() == 'cuda': + Tpetra.initialize_Kokkos(device_id=rank) + else: + Tpetra.initialize_Kokkos(num_threads=12) + # Use a seperate function to make sure that all Tpetra objects get deallocated before we finalize Kokkos. + main() + # finalize kokkos + Tpetra.finalize_Kokkos() diff --git a/packages/framework/ini-files/config-specs.ini b/packages/framework/ini-files/config-specs.ini index d7e4d008440a..e872329eac1e 100644 --- a/packages/framework/ini-files/config-specs.ini +++ b/packages/framework/ini-files/config-specs.ini @@ -781,6 +781,11 @@ opt-set-cmake-var CMAKE_CXX_COMPILER FILEPATH : ${SERIAL_CXX|ENV} opt-set-cmake-var CMAKE_Fortran_COMPILER FILEPATH : ${SERIAL_FC|ENV} opt-set-cmake-var TPL_ENABLE_MPI BOOL : OFF +[USE-MPI-NO-EXPLICIT-COMPILERS|YES] +opt-set-cmake-var TPL_ENABLE_MPI BOOL : ON + +[USE-MPI-NO-EXPLICIT-COMPILERS|NO] +opt-set-cmake-var TPL_ENABLE_MPI BOOL : OFF # @@ -1755,7 +1760,7 @@ use KOKKOS-ARCH|NO-KOKKOS-ARCH use USE-ASAN|NO use USE-FPIC|NO -use USE-MPI|YES +use USE-MPI-NO-EXPLICIT-COMPILERS|YES use USE-PT|NO use USE-COMPLEX|YES use USE-RDC|NO @@ -1805,7 +1810,7 @@ use KOKKOS-ARCH|NO-KOKKOS-ARCH use USE-ASAN|NO use USE-COMPLEX|NO use USE-FPIC|NO -use USE-MPI|NO +use USE-MPI-NO-EXPLICIT-COMPILERS|NO use USE-PT|NO use USE-RDC|NO use USE-UVM|NO @@ -1976,7 +1981,7 @@ opt-set-cmake-var TPL_FIND_SHARED_LIBS BOOL : OFF use KOKKOS-ARCH|AMPERE80 use USE-ASAN|NO use USE-FPIC|NO -use USE-MPI|YES +use USE-MPI-NO-EXPLICIT-COMPILERS|YES use USE-PT|NO use USE-COMPLEX|YES use USE-RDC|NO @@ -2003,7 +2008,7 @@ opt-set-cmake-var TPL_FIND_SHARED_LIBS BOOL : OFF use KOKKOS-ARCH|AMPERE80 use USE-ASAN|NO use USE-FPIC|NO -use USE-MPI|YES +use USE-MPI-NO-EXPLICIT-COMPILERS|YES use USE-PT|NO use USE-COMPLEX|YES use USE-RDC|NO diff --git a/packages/intrepid2/src/Projection/Intrepid2_ProjectionTools.hpp b/packages/intrepid2/src/Projection/Intrepid2_ProjectionTools.hpp index 79664a65f089..c7485e4a3563 100644 --- a/packages/intrepid2/src/Projection/Intrepid2_ProjectionTools.hpp +++ b/packages/intrepid2/src/Projection/Intrepid2_ProjectionTools.hpp @@ -50,7 +50,11 @@ #ifdef HAVE_INTREPID2_KOKKOSKERNELS #include "KokkosBatched_QR_Serial_Internal.hpp" #include "KokkosBatched_ApplyQ_Serial_Internal.hpp" +#if KOKKOS_VERSION >= 40599 +#include "KokkosBatched_Trsv_Decl.hpp" +#else #include "KokkosBatched_Trsv_Serial_Internal.hpp" +#endif #include "KokkosBatched_Util.hpp" #endif @@ -545,11 +549,15 @@ class ProjectionTools { w.data()); // R0^{-1} b -> b +#if KOKKOS_VERSION >= 40599 + KokkosBatched::SerialTrsv::invoke(1.0, A0, b); +#else KokkosBatched::SerialTrsvInternalUpper::invoke(false, A0.extent(0), 1.0, A0.data(), A0.stride_0(), A0.stride_1(), b.data(), b.stride_0()); +#endif //scattering b into the basis coefficients for(ordinal_type i=0; i b +#if KOKKOS_VERSION >= 40599 + KokkosBatched::SerialTrsv::invoke(1.0, A, b); +#else KokkosBatched::SerialTrsvInternalUpper::invoke(false, A.extent(0), 1.0, A.data(), A.stride_0(), A.stride_1(), b.data(), b.stride_0()); +#endif //scattering b into the basis coefficients for(ordinal_type i=0; i #include +#include namespace MueLu { @@ -244,3 +248,5 @@ class TpetraOperatorAsRowMatrix : public Tpetra::RowMatrix -#include -#include - -#include "MueLu_TestHelpers.hpp" -#include "MueLu_Version.hpp" - -#include "MueLu_MLParameterListInterpreter.hpp" -#include "MueLu_Exceptions.hpp" - -namespace MueLuTests { - -TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(MLParameterListInterpreter, SetParameterList, Scalar, LocalOrdinal, GlobalOrdinal, Node) { -#include "MueLu_UseShortNames.hpp" - MUELU_TESTING_SET_OSTREAM; - MUELU_TESTING_LIMIT_SCOPE(Scalar, GlobalOrdinal, NO); - if (!TYPE_EQUAL(SC, double)) { - out << "Skipping for SC != double" << std::endl; - return; - } - out << "version: " << MueLu::Version() << std::endl; - - // TODO: this test can be done at compilation time -#if !defined(HAVE_MUELU_EPETRA) or !defined(HAVE_MUELU_IFPACK) or !defined(HAVE_MUELU_AMESOS) - MUELU_TESTING_DO_NOT_TEST(Xpetra::UseEpetra, "Epetra, Ifpack, Amesos"); -#endif - -#if !defined(HAVE_MUELU_IFPACK2) or !defined(HAVE_MUELU_AMESOS2) - MUELU_TESTING_DO_NOT_TEST(Xpetra::UseTpetra, "Tpetra, Ifpack2, Amesos2"); -#endif - - RCP A = TestHelpers::TestFactory::Build1DPoisson(99); - Teuchos::ParameterList galeriParameters; - galeriParameters.set("nx", 99); - RCP coordinates = Galeri::Xpetra::Utils::CreateCartesianCoordinates("1D", A->getRowMap(), galeriParameters); - - ArrayRCP fileList = TestHelpers::GetFileList(std::string("ParameterList/MLParameterListInterpreter/"), std::string(".xml")); - - for (int i = 0; i < fileList.size(); i++) { - out << "Processing file: " << fileList[i] << std::endl; - Teuchos::ParameterList myList; - myList.set("xml parameter file", "ParameterList/MLParameterListInterpreter/" + fileList[i]); - - Teuchos::ArrayRCP xcoord = coordinates->getDataNonConst(0); - myList.set("x-coordinates", xcoord.get()); - - MLParameterListInterpreter mueluFactory(myList, A->getRowMap()->getComm()); - - RCP H = mueluFactory.CreateHierarchy(); - H->GetLevel(0)->template Set >("A", A); - - mueluFactory.SetupHierarchy(*H); - - // TODO: check no unused parameters - // TODO: check results of Iterate() - } -} - -#define MUELU_ETI_GROUP(SC, LO, GO, Node) \ - TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(MLParameterListInterpreter, SetParameterList, SC, LO, GO, Node) - -#include - -} // namespace MueLuTests - -// TODO: some tests of the parameter list parser can be done without building the Hierarchy. diff --git a/packages/muelu/test/unit_tests/ParameterList/MLParameterListInterpreter/CMakeLists.txt b/packages/muelu/test/unit_tests/ParameterList/MLParameterListInterpreter/CMakeLists.txt deleted file mode 100644 index ca8348a3c1ba..000000000000 --- a/packages/muelu/test/unit_tests/ParameterList/MLParameterListInterpreter/CMakeLists.txt +++ /dev/null @@ -1,10 +0,0 @@ -# Note about the use of wildcard in CMakeLists.txt: CMake dont't know -# when new files is added. You need to re-run CMake manually to -# regenerate a build system incorporating the new file. -# YOU MUST ALSO TOUCH A CMAKE CONFIGURATION FILE WHEN YOU PUSH THE NEW -# FILE TO FORCE THE RECONFIGURE ON OTHER PEOPLE'S BUILDS. -FILE(GLOB xmlFiles RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.xml) - -TRIBITS_COPY_FILES_TO_BINARY_DIR(ParameterList_MLParameterListInterpreter_cp - SOURCE_FILES ${xmlFiles} - ) diff --git a/packages/muelu/test/unit_tests/ParameterList/MLParameterListInterpreter/ml_SmootherParamAtTopLevel.xml b/packages/muelu/test/unit_tests/ParameterList/MLParameterListInterpreter/ml_SmootherParamAtTopLevel.xml deleted file mode 100644 index 8e53790faf4d..000000000000 --- a/packages/muelu/test/unit_tests/ParameterList/MLParameterListInterpreter/ml_SmootherParamAtTopLevel.xml +++ /dev/null @@ -1,5 +0,0 @@ - - - - - diff --git a/packages/muelu/test/unit_tests/ParameterList/MLParameterListInterpreter/ml_SmootherParamMixed.xml b/packages/muelu/test/unit_tests/ParameterList/MLParameterListInterpreter/ml_SmootherParamMixed.xml deleted file mode 100644 index 05388abb6a8d..000000000000 --- a/packages/muelu/test/unit_tests/ParameterList/MLParameterListInterpreter/ml_SmootherParamMixed.xml +++ /dev/null @@ -1,14 +0,0 @@ - - - - - - - - - - - - - - \ No newline at end of file diff --git a/packages/muelu/test/unit_tests/ParameterList/MLParameterListInterpreter/ml_SmootherParamPerLevel.xml b/packages/muelu/test/unit_tests/ParameterList/MLParameterListInterpreter/ml_SmootherParamPerLevel.xml deleted file mode 100644 index c2af11b7177b..000000000000 --- a/packages/muelu/test/unit_tests/ParameterList/MLParameterListInterpreter/ml_SmootherParamPerLevel.xml +++ /dev/null @@ -1,7 +0,0 @@ - - - - - - - \ No newline at end of file diff --git a/packages/muelu/test/unit_tests/ParameterList/MLParameterListInterpreter/ml_SmootherParamSubLists.xml b/packages/muelu/test/unit_tests/ParameterList/MLParameterListInterpreter/ml_SmootherParamSubLists.xml deleted file mode 100644 index dd4a963c5b76..000000000000 --- a/packages/muelu/test/unit_tests/ParameterList/MLParameterListInterpreter/ml_SmootherParamSubLists.xml +++ /dev/null @@ -1,11 +0,0 @@ - - - - - - - - - - - \ No newline at end of file diff --git a/packages/panzer/disc-fe/src/Panzer_AssemblyEngine_impl.hpp b/packages/panzer/disc-fe/src/Panzer_AssemblyEngine_impl.hpp index 2c0c506ec771..e4948dfa3c81 100644 --- a/packages/panzer/disc-fe/src/Panzer_AssemblyEngine_impl.hpp +++ b/packages/panzer/disc-fe/src/Panzer_AssemblyEngine_impl.hpp @@ -364,6 +364,8 @@ evaluateBCs(const panzer::BCType bc_type, workset.alpha = in.alpha; workset.beta = betaValue; workset.time = in.time; + workset.step_size = in.step_size; + workset.stage_number = in.stage_number; workset.gather_seeds = in.gather_seeds; workset.evaluate_transient_terms = in.evaluate_transient_terms; diff --git a/packages/panzer/disc-fe/src/Panzer_InitialCondition_Builder.cpp b/packages/panzer/disc-fe/src/Panzer_InitialCondition_Builder.cpp index 2636e1b246f6..887f7a84357e 100644 --- a/packages/panzer/disc-fe/src/Panzer_InitialCondition_Builder.cpp +++ b/packages/panzer/disc-fe/src/Panzer_InitialCondition_Builder.cpp @@ -151,7 +151,9 @@ evaluateInitialCondition(WorksetContainer & wkstContainer, const std::map< std::string,Teuchos::RCP< PHX::FieldManager > >& phx_ic_field_managers, Teuchos::RCP loc, const panzer::LinearObjFactory& lo_factory, - const double time_stamp) + const double time_stamp, + const double step_size, + const int stage_number) { typedef LinearObjContainer LOC; panzer::Traits::PED ped; @@ -185,6 +187,8 @@ evaluateInitialCondition(WorksetContainer & wkstContainer, for (std::size_t i = 0; i < w.size(); ++i) { panzer::Workset& workset = w[i]; workset.time = time_stamp; + workset.step_size = step_size; + workset.stage_number = static_cast(stage_number); fm->evaluateFields(workset); } diff --git a/packages/panzer/disc-fe/src/Panzer_InitialCondition_Builder.hpp b/packages/panzer/disc-fe/src/Panzer_InitialCondition_Builder.hpp index a0b4ddb71e1b..551ccbbe9533 100644 --- a/packages/panzer/disc-fe/src/Panzer_InitialCondition_Builder.hpp +++ b/packages/panzer/disc-fe/src/Panzer_InitialCondition_Builder.hpp @@ -117,12 +117,16 @@ namespace panzer { std::vector > & physics_blocks); /** Execute the construction of an initial condition. + * + * The time_stamp, step_size and stage_number are workset parameters used in some transient problems. These are optional. */ void evaluateInitialCondition(WorksetContainer & wkstContainer, const std::map > >& phx_ic_field_managers, Teuchos::RCP loc, const panzer::LinearObjFactory& lo_factory, - const double time_stamp); + const double time_stamp, + const double step_size = 0.0, + const int stage_number = 0); } #endif diff --git a/packages/piro/src/Piro_ObserverToTempusIntegrationObserverAdapter.hpp b/packages/piro/src/Piro_ObserverToTempusIntegrationObserverAdapter.hpp index 6c23c783c6ba..4a4bccea08d1 100644 --- a/packages/piro/src/Piro_ObserverToTempusIntegrationObserverAdapter.hpp +++ b/packages/piro/src/Piro_ObserverToTempusIntegrationObserverAdapter.hpp @@ -40,31 +40,13 @@ class ObserverToTempusIntegrationObserverAdapter : public Tempus::IntegratorObse //@{ /// Destructor - virtual ~ObserverToTempusIntegrationObserverAdapter(); + virtual ~ObserverToTempusIntegrationObserverAdapter() = default; /// Observe the beginning of the time integrator. virtual void observeStartIntegrator(const Tempus::Integrator& integrator) override; - /// Observe the beginning of the time step loop. - virtual void observeStartTimeStep(const Tempus::Integrator& integrator) override; - - /// Observe after the next time step size is selected. - virtual void observeNextTimeStep(const Tempus::Integrator& integrator) override; - - /// Observe before Stepper takes step. - virtual void observeBeforeTakeStep(const Tempus::Integrator& integrator) override; - - /// Observe after Stepper takes step. - virtual void observeAfterTakeStep(const Tempus::Integrator& integrator) override; - - /// Observe after checking time step. Observer can still fail the time step here. - virtual void observeAfterCheckTimeStep(const Tempus::Integrator& integrator) override; - /// Observe the end of the time step loop. virtual void observeEndTimeStep(const Tempus::Integrator& integrator) override; - - /// Observe the end of the time integrator. - virtual void observeEndIntegrator(const Tempus::Integrator& integrator) override; //@} private: diff --git a/packages/piro/src/Piro_ObserverToTempusIntegrationObserverAdapter_Def.hpp b/packages/piro/src/Piro_ObserverToTempusIntegrationObserverAdapter_Def.hpp index 441d815b13bd..8724eabe6275 100644 --- a/packages/piro/src/Piro_ObserverToTempusIntegrationObserverAdapter_Def.hpp +++ b/packages/piro/src/Piro_ObserverToTempusIntegrationObserverAdapter_Def.hpp @@ -32,12 +32,6 @@ Piro::ObserverToTempusIntegrationObserverAdapter::ObserverToTempusIntegr previous_dt_ = 0.0; } -template -Piro::ObserverToTempusIntegrationObserverAdapter::~ObserverToTempusIntegrationObserverAdapter() -{ - //Nothing to do -} - template void Piro::ObserverToTempusIntegrationObserverAdapter:: @@ -64,49 +58,6 @@ observeStartIntegrator(const Tempus::Integrator& integrator) this->observeTimeStep(); } -template -void -Piro::ObserverToTempusIntegrationObserverAdapter:: -observeStartTimeStep(const Tempus::Integrator& ) -{ - //Nothing to do -} - -template -void -Piro::ObserverToTempusIntegrationObserverAdapter:: -observeNextTimeStep(const Tempus::Integrator& ) -{ - //Nothing to do -} - -template -void -Piro::ObserverToTempusIntegrationObserverAdapter:: -observeBeforeTakeStep(const Tempus::Integrator& ) -{ - //Nothing to do -} - - -template -void -Piro::ObserverToTempusIntegrationObserverAdapter:: -observeAfterTakeStep(const Tempus::Integrator& ) -{ - //Nothing to do -} - - -template -void -Piro::ObserverToTempusIntegrationObserverAdapter:: -observeAfterCheckTimeStep(const Tempus::Integrator& integrator) -{ - //Nothing to do -} - - template void Piro::ObserverToTempusIntegrationObserverAdapter:: @@ -139,35 +90,6 @@ observeEndTimeStep(const Tempus::Integrator& integrator) this->observeTimeStep(); } - -template -void -Piro::ObserverToTempusIntegrationObserverAdapter:: -observeEndIntegrator(const Tempus::Integrator& integrator) -{ - //this->observeTimeStep(); - - std::string exitStatus; - //const Scalar runtime = integrator.getIntegratorTimer()->totalElapsedTime(); - if (integrator.getSolutionHistory()->getCurrentState()->getSolutionStatus() == - Tempus::Status::FAILED or integrator.getStatus() == Tempus::Status::FAILED) { - exitStatus = "Time integration FAILURE!"; - } else { - exitStatus = "Time integration complete."; - } - std::time_t end = std::time(nullptr); - const Scalar runtime = integrator.getIntegratorTimer()->totalElapsedTime(); - const Teuchos::RCP out = integrator.getOStream(); - Teuchos::OSTab ostab(out,0,"ScreenOutput"); - *out << "============================================================================\n" - << " Total runtime = " << runtime << " sec = " - << runtime/60.0 << " min\n" - << std::asctime(std::localtime(&end)) - << exitStatus << "\n" - << std::endl; - -} - template void Piro::ObserverToTempusIntegrationObserverAdapter::observeTimeStep() diff --git a/packages/piro/src/Piro_TempusIntegrator.hpp b/packages/piro/src/Piro_TempusIntegrator.hpp index 38e4928e90a4..07b677f7b8e2 100644 --- a/packages/piro/src/Piro_TempusIntegrator.hpp +++ b/packages/piro/src/Piro_TempusIntegrator.hpp @@ -82,6 +82,7 @@ class TempusIntegrator //The following routine is only for adjoint sensitivities Teuchos::RCP> getDgDp() const; + Teuchos::RCP > getIntegrator () const; private: Teuchos::RCP > basicIntegrator_; diff --git a/packages/piro/src/Piro_TempusIntegrator_Def.hpp b/packages/piro/src/Piro_TempusIntegrator_Def.hpp index beb0d58be185..fb25f7ddf80d 100644 --- a/packages/piro/src/Piro_TempusIntegrator_Def.hpp +++ b/packages/piro/src/Piro_TempusIntegrator_Def.hpp @@ -363,3 +363,22 @@ Piro::TempusIntegrator::getDgDp() const "which is not of type Tempus::IntegratorAdjointSensitivity!\n"); } } + + +template +Teuchos::RCP > +Piro::TempusIntegrator:: getIntegrator () const +{ + Teuchos::RCP > out; + if (basicIntegrator_ != Teuchos::null) { + out = basicIntegrator_; + } else if (fwdSensIntegrator_ != Teuchos::null) { + out = fwdSensIntegrator_; + } else if (adjSensIntegrator_ != Teuchos::null) { + out = adjSensIntegrator_; + } else { + TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, + "Error in Piro::TempusIntegrator: no integrator stored!\n"); + } + return out; +} diff --git a/packages/piro/src/Piro_TempusSolver.hpp b/packages/piro/src/Piro_TempusSolver.hpp index 6912d13806cd..fd31ba53f76a 100644 --- a/packages/piro/src/Piro_TempusSolver.hpp +++ b/packages/piro/src/Piro_TempusSolver.hpp @@ -161,7 +161,6 @@ class TempusSolver Teuchos::RCP > model_; Teuchos::RCP > adjointModel_; - Teuchos::RCP > thyraModel_; Scalar t_initial_; Scalar t_final_; diff --git a/packages/piro/src/Piro_TempusSolver_Def.hpp b/packages/piro/src/Piro_TempusSolver_Def.hpp index 333e9824705b..e11a298d200e 100644 --- a/packages/piro/src/Piro_TempusSolver_Def.hpp +++ b/packages/piro/src/Piro_TempusSolver_Def.hpp @@ -167,24 +167,10 @@ void Piro::TempusSolver::initialize( const std::string stepperType = stepperPL->get("Stepper Type", "Backward Euler"); //*out_ << "Stepper Type = " << stepperType << "\n"; - // - // *out_ << "\nB) Create the Stratimikos linear solver factory ...\n"; - // - // This is the linear solve strategy that will be used to solve for the - // linear system with the W. - // - Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder; - -#ifdef HAVE_PIRO_MUELU - Stratimikos::enableMueLu(linearSolverBuilder); -#endif - - linearSolverBuilder.setParameterList(sublist(tempusPL, "Stratimikos", true)); tempusPL->validateParameters(*getValidTempusParameters( tempusPL->get("Integrator Name", "Tempus Integrator"), integratorPL->get("Stepper Name", "Tempus Stepper") ), 0); - RCP > lowsFactory = createLinearSolveStrategy(linearSolverBuilder); // *out_ << "\nC) Create and initalize the forward model ...\n"; @@ -287,11 +273,6 @@ void Piro::TempusSolver::initialize( } } - // C.2) Create the Thyra-wrapped ModelEvaluator - - thyraModel_ = rcp(new Thyra::DefaultModelEvaluatorWithSolveFactory(model_, lowsFactory)); - const RCP > x_space = thyraModel_->get_x_space(); - // *out_ << "\nD) Create the stepper and integrator for the forward problem ...\n"; @@ -650,8 +631,12 @@ template void Piro::TempusSolver:: setObserver() const { - Teuchos::RCP > observer = Teuchos::null; - if (Teuchos::nonnull(piroObserver_)) { + // Do not create the tempus observer adapter if the user-provided Piro observer is already a tempus observer + Teuchos::RCP > observer; + observer = Teuchos::rcp_dynamic_cast>(piroObserver_); + if (Teuchos::is_null(observer) and Teuchos::nonnull(piroObserver_)) { + // The user did not provide a Tempus observer, so create an adapter one + //Get solutionHistory from integrator const Teuchos::RCP > solutionHistory = piroTempusIntegrator_->getSolutionHistory(); const Teuchos::RCP > timeStepControl = piroTempusIntegrator_->getTimeStepControl(); @@ -659,6 +644,7 @@ setObserver() const observer = Teuchos::rcp(new ObserverToTempusIntegrationObserverAdapter(solutionHistory, timeStepControl, piroObserver_, supports_x_dotdot_, abort_on_fail_at_min_dt_, sens_method_)); } + if (Teuchos::nonnull(observer)) { //Set observer in integrator piroTempusIntegrator_->clearObservers(); diff --git a/packages/piro/src/Piro_TransientSolver_Def.hpp b/packages/piro/src/Piro_TransientSolver_Def.hpp index 358bae08365e..bd303be53a14 100644 --- a/packages/piro/src/Piro_TransientSolver_Def.hpp +++ b/packages/piro/src/Piro_TransientSolver_Def.hpp @@ -348,7 +348,17 @@ Piro::TransientSolver::evalConvergedModelResponsesAndSensitivities( // Solution at convergence is the response at index num_g_ RCP > gx_out = outArgs.get_g(num_g_); if (Teuchos::nonnull(gx_out)) { - Thyra::copy(*modelInArgs.get_x(), gx_out.ptr()); + auto x = modelInArgs.get_x(); + if (x->space()->isCompatible(*gx_out->space())) { + Thyra::copy(*modelInArgs.get_x(), gx_out.ptr()); + } else { + *out_ << " WARNING [PIRO::TransientSolver::evalConvergedModelResponsesAndSensitivities]\n" + " The solution from inArgs (x) is incompatible with the last response in outArgs (gx),\n" + " which was created as a vector with the same vector-space of x. Since the responses\n" + " were created BEFORE calling the transient solver, the most likely explanation for this\n" + " incompatibility is that during the time integration the solution vector space changed,\n" + " for instance because the underlying spatial mesh was adapted (with topological changes).\n"; + } } // Setup output for final evalution of underlying model diff --git a/packages/piro/test/CMakeLists.txt b/packages/piro/test/CMakeLists.txt index b964bf3d61f7..0471915e9e9e 100644 --- a/packages/piro/test/CMakeLists.txt +++ b/packages/piro/test/CMakeLists.txt @@ -152,6 +152,19 @@ IF (Piro_ENABLE_NOX) ENDIF(Piro_ENABLE_NOX) IF (Piro_ENABLE_Tempus) + TRIBITS_ADD_EXECUTABLE_AND_TEST( + TempusSolver_UnitTests_Tpetra + SOURCES + Piro_TempusSolver_UnitTests.cpp + ${TEUCHOS_STD_UNIT_TEST_MAIN} + Piro_Test_ThyraSupport.hpp + MockModelEval_A_Tpetra.hpp + MockModelEval_A_Tpetra.cpp + MatrixBased_LOWS.cpp + NUM_MPI_PROCS 1-4 + STANDARD_PASS_OUTPUT + ) + TRIBITS_ADD_EXECUTABLE_AND_TEST( TempusSolver_SensitivitySinCos_Combined_FSA_UnitTests SOURCES @@ -251,6 +264,7 @@ TRIBITS_COPY_FILES_TO_BINARY_DIR(copyTestInputFiles input_Analysis_ROL_ReducedSpace_TrustRegion_BoundConstrained_ExplicitAdjointME_NOXSolver.xml input_Analysis_ROL_FullSpace_AugmentedLagrangian_BoundConstrained.xml input_Tempus_BackwardEuler_SinCos.xml + input_tempus_be_nox_solver.yaml SOURCE_DIR ${PACKAGE_SOURCE_DIR}/test SOURCE_PREFIX "_" EXEDEPS ${ThyraSolverTpetra_EXENAME} ${AnalysisDriverTpetra_EXENAME} @@ -380,6 +394,7 @@ IF (PIRO_HAVE_EPETRA_STACK) MockModelEval_A.cpp NUM_MPI_PROCS 1-4 STANDARD_PASS_OUTPUT + TARGET_DEFINES TEST_USE_EPETRA ) TRIBITS_ADD_EXECUTABLE_AND_TEST( TempusSolverForwardOnly_UnitTests @@ -408,7 +423,6 @@ IF (PIRO_HAVE_EPETRA_STACK) NUM_MPI_PROCS 4 STANDARD_PASS_OUTPUT ) - ENDIF (Piro_ENABLE_Tempus) TRIBITS_COPY_FILES_TO_BINARY_DIR(copyEpetraTestInputFiles diff --git a/packages/piro/test/MockModelEval_A_Tpetra.cpp b/packages/piro/test/MockModelEval_A_Tpetra.cpp index 54fcf701db58..a7cb776fe1ad 100644 --- a/packages/piro/test/MockModelEval_A_Tpetra.cpp +++ b/packages/piro/test/MockModelEval_A_Tpetra.cpp @@ -231,6 +231,7 @@ MockModelEval_A_Tpetra::createOutArgsImpl() const result.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_f, true); result.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_W_op, true); + result.setSupports(Thyra::ModelEvaluatorBase::OUT_ARG_W, true); result.set_W_properties(Thyra::ModelEvaluatorBase::DerivativeProperties( Thyra::ModelEvaluatorBase::DERIV_LINEARITY_UNKNOWN, Thyra::ModelEvaluatorBase::DERIV_RANK_FULL, diff --git a/packages/piro/test/Piro_TempusSolver_UnitTests.cpp b/packages/piro/test/Piro_TempusSolver_UnitTests.cpp index 16b18ffc952e..721fb4351b57 100644 --- a/packages/piro/test/Piro_TempusSolver_UnitTests.cpp +++ b/packages/piro/test/Piro_TempusSolver_UnitTests.cpp @@ -13,6 +13,7 @@ #include "Piro_TempusSolver.hpp" #include "Tempus_StepperFactory.hpp" #include "Piro_ObserverToTempusIntegrationObserverAdapter.hpp" +#include "Piro_Test_MockTempusObserver.hpp" #ifdef HAVE_PIRO_NOX #include "Piro_NOXSolver.hpp" @@ -23,15 +24,21 @@ #include "Piro_Test_MockObserver.hpp" #include "Piro_TempusIntegrator.hpp" -#include "MockModelEval_A.hpp" +#ifdef TEST_USE_EPETRA #include "Thyra_EpetraModelEvaluator.hpp" +#include "MockModelEval_A.hpp" +#include "Thyra_AmesosLinearOpWithSolveFactory.hpp" +#else +#include "MockModelEval_A_Tpetra.hpp" +#include "Piro_StratimikosUtils.hpp" +#include "Teuchos_YamlParameterListCoreHelpers.hpp" +#endif + #include "Thyra_ModelEvaluatorHelpers.hpp" #include "Thyra_DefaultNominalBoundsOverrideModelEvaluator.hpp" -#include "Thyra_AmesosLinearOpWithSolveFactory.hpp" - #include "Teuchos_UnitTestHarness.hpp" #include "Teuchos_Ptr.hpp" @@ -51,6 +58,7 @@ namespace Thyra { // Setup support +#ifdef TEST_USE_EPETRA const RCP epetraModelNew() { #ifdef HAVE_MPI @@ -60,20 +68,39 @@ const RCP epetraModelNew() #endif /*HAVE_MPI*/ return rcp(new MockModelEval_A(comm)); } - -const RCP > thyraModelNew(const RCP &epetraModel) +const RCP > thyraModelNew(const RCP &epetraModel) { const RCP > lowsFactory(new Thyra::AmesosLinearOpWithSolveFactory); return epetraModelEvaluator(epetraModel, lowsFactory); } -RCP > defaultModelNew() +RCP > defaultModelNew() { return thyraModelNew(epetraModelNew()); } +#else +RCP > defaultModelNew() +{ + auto comm = Tpetra::getDefaultComm(); + auto model = rcp(new MockModelEval_A_Tpetra(comm)); + + std::string inputFile = "input_tempus_be_nox_solver.yaml"; + auto piroParams = rcp(new Teuchos::ParameterList("Piro Parameters")); + Teuchos::updateParametersFromYamlFile(inputFile, piroParams.ptr()); + + auto stratParams = Piro::extractStratimikosParams(piroParams); + + Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder; + linearSolverBuilder.setParameterList(stratParams); + + auto lowsFactory = createLinearSolveStrategy(linearSolverBuilder); + + return rcp(new Thyra::DefaultModelEvaluatorWithSolveFactory(model, lowsFactory)); +} +#endif const RCP > solverNew( - const RCP > &thyraModel, + const RCP > &thyraModel, double finalTime, const std::string sens_method_string) { @@ -105,7 +132,7 @@ const RCP > solverNew( } const RCP > solverNew( - const RCP > &thyraModel, + const RCP > &thyraModel, double initialTime, double finalTime, const RCP > &observer, @@ -145,7 +172,7 @@ const RCP > solverNew( } const RCP > solverNew( - const RCP > &thyraModel, + const RCP > &thyraModel, double initialTime, double finalTime, double fixedTimeStep, @@ -167,6 +194,7 @@ const RCP > solverNew( tempusPL->sublist("Demo Stepper").set("Zero Initial Guess", false); tempusPL->sublist("Demo Stepper").set("Solver Name", "Demo Solver"); tempusPL->sublist("Demo Stepper").sublist("Demo Solver").sublist("NOX").sublist("Direction").set("Method","Newton"); + Teuchos::RCP > integrator = Teuchos::rcp(new Piro::TempusIntegrator(tempusPL, thyraModel)); const RCP > solutionHistory = integrator->getSolutionHistory(); const RCP > timeStepControl = integrator->getTimeStepControl(); @@ -202,7 +230,7 @@ const double tol = 1.0e-8; TEUCHOS_UNIT_TEST(Piro_TempusSolver, TimeZero_Solution) { - const RCP > model = defaultModelNew(); + const auto model = defaultModelNew(); const double finalTime = 0.0; const RCP > solver = solverNew(model, finalTime, "None"); @@ -229,7 +257,7 @@ TEUCHOS_UNIT_TEST(Piro_TempusSolver, TimeZero_Solution) TEUCHOS_UNIT_TEST(Piro_TempusSolver, TimeZero_Response) { - const RCP > model = defaultModelNew(); + const auto model = defaultModelNew(); const int responseIndex = 0; @@ -305,7 +333,7 @@ TEUCHOS_UNIT_TEST(Piro_TempusSolver, TimeZero_NoDgDp_NoResponseSensitivity) TEUCHOS_UNIT_TEST(Piro_TempusSolver, ObserveInitialCondition) { - const RCP > model = defaultModelNew(); + const auto model = defaultModelNew(); const RCP > observer(new MockObserver); const double timeStamp = 2.0; @@ -334,15 +362,100 @@ TEUCHOS_UNIT_TEST(Piro_TempusSolver, ObserveInitialCondition) TEUCHOS_UNIT_TEST(Piro_TempusSolver, ObserveFinalSolution) { - const RCP > model = defaultModelNew(); + const auto model = defaultModelNew(); const RCP > observer(new MockObserver); const double initialTime = 0.0; const double finalTime = 0.1; const double timeStepSize = 0.05; - const RCP > solver = - solverNew(model, initialTime, finalTime, timeStepSize, observer, "None"); + const RCP appParams(new ParameterList("params")); +#ifdef TEST_USE_EPETRA + auto& tempusPL = appParams->sublist("Tempus"); + tempusPL.set("Integrator Name", "Demo Integrator"); + tempusPL.sublist("Demo Integrator").set("Integrator Type", "Integrator Basic"); + tempusPL.sublist("Demo Integrator").set("Stepper Name", "Demo Stepper"); + tempusPL.sublist("Demo Integrator").sublist("Solution History").set("Storage Type", "Unlimited"); + tempusPL.sublist("Demo Integrator").sublist("Solution History").set("Storage Limit", 20); + tempusPL.sublist("Demo Stepper").set("Stepper Type", "Backward Euler"); + tempusPL.sublist("Demo Stepper").set("Zero Initial Guess", false); + tempusPL.sublist("Demo Stepper").set("Solver Name", "Demo Solver"); + tempusPL.sublist("Demo Stepper").sublist("Demo Solver").sublist("NOX").sublist("Direction").set("Method","Newton"); +#else + Teuchos::updateParametersFromYamlFile("input_tempus_be_nox_solver.yaml",appParams.ptr()); +#endif + auto& tsc = appParams->sublist("Tempus").sublist("Demo Integrator").sublist("Time Step Control"); + tsc.set("Initial Time", initialTime); + tsc.set("Final Time", finalTime); + tsc.set("Minimum Time Step", timeStepSize/10); + tsc.set("Initial Time Step", timeStepSize); + tsc.set("Maximum Time Step", timeStepSize); + + RCP> adjointModel = Teuchos::null; + RCP> piro_obs = observer; + auto solver = tempusSolver(appParams,model,adjointModel,piro_obs); + + const Thyra::MEB::InArgs inArgs = solver->getNominalValues(); + + Thyra::MEB::OutArgs outArgs = solver->createOutArgs(); + const int solutionResponseIndex = solver->Ng() - 1; + const RCP > solution = + Thyra::createMember(solver->get_g_space(solutionResponseIndex)); + outArgs.set_g(solutionResponseIndex, solution); + + solver->evalModel(inArgs, outArgs); + + TEST_COMPARE_FLOATING_ARRAYS( + arrayFromVector(*observer->lastSolution()), + arrayFromVector(*solution), + tol); + + TEST_FLOATING_EQUALITY(observer->lastStamp(), finalTime, tol); +} + +TEUCHOS_UNIT_TEST(Piro_TempusSolver, ObserveFinalSolutionNoWrapper) +{ + const auto model = defaultModelNew(); + auto observer = Teuchos::rcp(new MockTempusObserver()); + + const double initialTime = 0.0; + const double finalTime = 0.1; + const double timeStepSize = 0.05; + const RCP appParams(new ParameterList("params")); +#ifdef TEST_USE_EPETRA + auto& tempusPL = appParams->sublist("Tempus"); + tempusPL.set("Integrator Name", "Demo Integrator"); + tempusPL.sublist("Demo Integrator").set("Integrator Type", "Integrator Basic"); + tempusPL.sublist("Demo Integrator").set("Stepper Name", "Demo Stepper"); + tempusPL.sublist("Demo Integrator").sublist("Solution History").set("Storage Type", "Unlimited"); + tempusPL.sublist("Demo Integrator").sublist("Solution History").set("Storage Limit", 20); + tempusPL.sublist("Demo Stepper").set("Stepper Type", "Backward Euler"); + tempusPL.sublist("Demo Stepper").set("Zero Initial Guess", false); + tempusPL.sublist("Demo Stepper").set("Solver Name", "Demo Solver"); + tempusPL.sublist("Demo Stepper").sublist("Demo Solver").sublist("NOX").sublist("Direction").set("Method","Newton"); +#else + Teuchos::updateParametersFromYamlFile("input_tempus_be_nox_solver.yaml",appParams.ptr()); +#endif + auto& tsc = appParams->sublist("Tempus").sublist("Demo Integrator").sublist("Time Step Control"); + tsc.set("Initial Time", initialTime); + tsc.set("Final Time", finalTime); + tsc.set("Minimum Time Step", timeStepSize/10); + tsc.set("Initial Time Step", timeStepSize); + tsc.set("Maximum Time Step", timeStepSize); + + RCP> adjointModel = Teuchos::null; + RCP> piro_obs = observer; + auto solver = tempusSolver(appParams,model,adjointModel,piro_obs); + + // Check that the piro integrator did NOT wrap our observer in a tempus adapter + auto piro_integrator = solver->getPiroTempusIntegrator(); + auto tempus_integrator = piro_integrator->getIntegrator(); + auto basic_integrator = Teuchos::rcp_dynamic_cast>(tempus_integrator); + auto nonconst_basic_integrator = Teuchos::rcp_const_cast>(basic_integrator); + auto tempus_observer = nonconst_basic_integrator->getObserver(); + TEST_ASSERT (tempus_observer==observer); + + // Proceed like the test above const Thyra::MEB::InArgs inArgs = solver->getNominalValues(); Thyra::MEB::OutArgs outArgs = solver->createOutArgs(); diff --git a/packages/piro/test/Piro_Test_MockTempusObserver.hpp b/packages/piro/test/Piro_Test_MockTempusObserver.hpp new file mode 100644 index 000000000000..40f528439ca0 --- /dev/null +++ b/packages/piro/test/Piro_Test_MockTempusObserver.hpp @@ -0,0 +1,68 @@ +// @HEADER +// ***************************************************************************** +// Piro: Strategy package for embedded analysis capabilitites +// +// Copyright 2010 NTESS and the Piro contributors. +// SPDX-License-Identifier: BSD-3-Clause +// ***************************************************************************** +// @HEADER + +#ifndef PIRO_TEST_MOCKTEMPUSOBSERVER_HPP +#define PIRO_TEST_MOCKTEMPUSOBSERVER_HPP + +#include "Piro_Test_MockObserver.hpp" +#include "Tempus_IntegratorObserverBasic.hpp" + +namespace Piro { + +namespace Test { + +template +class MockTempusObserver : public MockObserver, + public Tempus::IntegratorObserverBasic +{ +public: + MockTempusObserver() = default; + + void observeEndTimeStep(const Tempus::Integrator& integrator) override; +}; + +// ------------ IMPL ------------ // + +template +void MockTempusObserver:: +observeEndTimeStep(const Tempus::Integrator& integrator) +{ + // Simply grab stuff from integrator and forward to the base class observeSolution method + // NOTE: this simple class IGNORES completely dxdp + + auto sol_hist = integrator.getSolutionHistory(); + + auto x = sol_hist->getCurrentState()->getX(); + auto x_dot = sol_hist->getCurrentState()->getXDot(); + auto x_dotdot = sol_hist->getCurrentState()->getXDotDot(); + + const Scalar scalar_time = sol_hist->getCurrentState()->getTime(); + typedef typename Teuchos::ScalarTraits::magnitudeType StampScalar; + const StampScalar time = Teuchos::ScalarTraits::real(scalar_time); + + if (x_dot != Teuchos::null) { + if (x_dotdot != Teuchos::null) { + // x_dot AND x_dotdot + this->observeSolution(*x, *x_dot, *x_dotdot, time); + } else { + // no x_dotdot + this->observeSolution(*x, *x_dot, time); + } + } else { + //no x_dot + this->observeSolution(*x, time); + } +} + +} // namespace Test + +} // namespace Piro + +#endif /* PIRO_TEST_MOCKTEMPUSOBSERVER_HPP */ + diff --git a/packages/piro/test/_input_tempus_be_nox_solver.yaml b/packages/piro/test/_input_tempus_be_nox_solver.yaml new file mode 100644 index 000000000000..9f0d8dd35231 --- /dev/null +++ b/packages/piro/test/_input_tempus_be_nox_solver.yaml @@ -0,0 +1,90 @@ +Piro: + Solver Type: Tempus + Tempus: + Integrator Name: Demo Integrator + Demo Integrator: + Integrator Type: Integrator Basic + Screen Output Index List: '1' + Screen Output Index Interval: 100 + Stepper Name: Demo Stepper + Solution History: + Storage Type: Unlimited + Storage Limit: 20 + Time Step Control: + Initial Time: 0.0 + Initial Time Index: 0 + Initial Time Step: 5.0e-03 + Final Time: 8.0e-1 + Final Time Index: 10000 + Maximum Absolute Error: 1.0e-08 + Maximum Relative Error: 1.0e-08 + Output Time List: '' + Output Index List: '' + Output Index Interval: 1000 + Maximum Number of Stepper Failures: 10 + Maximum Number of Consecutive Stepper Failures: 5 + Demo Stepper: + Stepper Type: Backward Euler + Solver Name: Demo Solver + Demo Solver: + NOX: + Direction: + Method: Newton + Newton: + Forcing Term Method: Constant + Rescue Bad Newton Solve: true + Line Search: + Full Step: + Full Step: 1.0e+00 + Method: Backtrack + Nonlinear Solver: Line Search Based + Printing: + Output Precision: 3 + Output Processor: 0 + Output Information: + Error: true + Warning: true + Outer Iteration: false + Parameters: true + Details: false + Linear Solver Details: true + Stepper Iteration: true + Stepper Details: true + Stepper Parameters: true + Solver Options: + Status Test Check Type: Minimal + Status Tests: + Test Type: Combo + Combo Type: OR + Number of Tests: 2 + Test 0: + Test Type: NormF + Tolerance: 1.0e-08 + Test 1: + Test Type: MaxIters + Maximum Iterations: 20 + Stratimikos: + Linear Solver Type: Belos + Linear Solver Types: + Belos: + Solver Type: Block GMRES + Solver Types: + Block GMRES: + Convergence Tolerance: 1.0e-05 + Output Frequency: 10 + Output Style: 1 + Verbosity: 33 + Maximum Iterations: 200 + Block Size: 1 + Num Blocks: 4 + Flexible Gmres: false + Preconditioner Type: Ifpack2 + Preconditioner Types: + Ifpack2: + Prec Type: RILUK + Overlap: 1 + Ifpack2 Settings: + 'fact: drop tolerance': 0.0 + 'fact: iluk level-of-fill': 0 + 'fact: level-of-fill': 1 +... diff --git a/packages/shylu/shylu_dd/frosch/src/InterfaceSets/FROSch_InterfaceEntity_def.hpp b/packages/shylu/shylu_dd/frosch/src/InterfaceSets/FROSch_InterfaceEntity_def.hpp index 1b842a939a75..8ec3e3b1f346 100644 --- a/packages/shylu/shylu_dd/frosch/src/InterfaceSets/FROSch_InterfaceEntity_def.hpp +++ b/packages/shylu/shylu_dd/frosch/src/InterfaceSets/FROSch_InterfaceEntity_def.hpp @@ -117,7 +117,7 @@ namespace FROSch { FROSCH_ASSERT(nDofs<=DofsPerNode_,"nDofs>DofsPerNode_."); for (unsigned i=0; i { /** \brief Set the number of time events. * - * - If numEvents_ < 0, then reset numEvents from start_, stop_ and stride_. - * - ElseIf start_ = stop_, numEvents = 1, and reset stride = 0. - * - ElseIf numEvents_ < 2, numEvents = 2, and stride = stop_ - start_. + * - If numEvents < 0, then reset numEvents_ from start_, stop_ and stride_. + * - ElseIf start_ = stop_, numEvents_ = 1, and reset stride = 0. + * - ElseIf numEvents < 2, numEvents_ = 2, and stride = stop_ - start_. * - Else stride = (stop_ - start_)/(numEvents_-1). * * If the resulting stride is less than twice the absolute tolerance, diff --git a/packages/tempus/src/Tempus_TimeEventRange_impl.hpp b/packages/tempus/src/Tempus_TimeEventRange_impl.hpp index 0781798154a5..867dfdf73aaa 100644 --- a/packages/tempus/src/Tempus_TimeEventRange_impl.hpp +++ b/packages/tempus/src/Tempus_TimeEventRange_impl.hpp @@ -169,19 +169,18 @@ void TimeEventRange::setTimeStride(Scalar stride) template void TimeEventRange::setNumEvents(int numEvents) { - numEvents_ = numEvents; - if (numEvents_ < 0) { // Do not use numEvents_ to set stride! Reset - // numEvents_ from stride_. + if (numEvents < 0) { // Do not use numEvents to set stride! + // Reset numEvents_ from stride_. if (stride_ < 2 * absTol_) stride_ = 2 * absTol_; numEvents_ = int((stop_ + absTol_ - start_) / stride_) + 1; return; } else if (approxEqualAbsTol(start_, stop_, absTol_)) { numEvents_ = 1; - stride_ = 0.0; stride_ = stop_ - start_; } else { + numEvents_ = numEvents; if (numEvents_ < 2) numEvents_ = 2; stride_ = (stop_ - start_) / Scalar(numEvents_ - 1); } diff --git a/packages/tempus/unit_test/Tempus_UnitTest_TimeEventRange.cpp b/packages/tempus/unit_test/Tempus_UnitTest_TimeEventRange.cpp index d4d978f1b54f..8bb0863edc28 100644 --- a/packages/tempus/unit_test/Tempus_UnitTest_TimeEventRange.cpp +++ b/packages/tempus/unit_test/Tempus_UnitTest_TimeEventRange.cpp @@ -154,12 +154,34 @@ TEUCHOS_UNIT_TEST(TimeEventRange, setTimeRange) { auto te = rcp(new Tempus::TimeEventRange()); - // Set with time range. + // Set with time range with stride. te->setTimeRange(0.0, PI, 1.0); TEST_FLOATING_EQUALITY(te->getTimeStart(), 0.0, 1.0e-14); TEST_FLOATING_EQUALITY(te->getTimeStop(), PI, 1.0e-14); TEST_FLOATING_EQUALITY(te->getTimeStride(), 1.0, 1.0e-14); TEST_COMPARE(te->getNumEvents(), ==, 4); + + // Set with time range number of events. + te->setTimeRange(0.0, PI, 5); + TEST_FLOATING_EQUALITY(te->getTimeStart(), 0.0, 1.0e-14); + TEST_FLOATING_EQUALITY(te->getTimeStop(), PI, 1.0e-14); + TEST_FLOATING_EQUALITY(te->getTimeStride(), PI / 4.0, 1.0e-14); + TEST_COMPARE(te->getNumEvents(), ==, 5); +} + +// ************************************************************ +// ************************************************************ +TEUCHOS_UNIT_TEST(TimeEventRange, setNumEvents) +{ + auto te = rcp(new Tempus::TimeEventRange()); + + // Set with time range with stride, then set number of events. + te->setTimeRange(0.0, PI, 1.0); + te->setNumEvents(5); + TEST_FLOATING_EQUALITY(te->getTimeStart(), 0.0, 1.0e-14); + TEST_FLOATING_EQUALITY(te->getTimeStop(), PI, 1.0e-14); + TEST_FLOATING_EQUALITY(te->getTimeStride(), PI / 4.0, 1.0e-14); + TEST_COMPARE(te->getNumEvents(), ==, 5); } // ************************************************************ diff --git a/packages/thyra/core/src/interfaces/operator_solve/fundamental/Thyra_LinearOpWithSolveBase_decl.hpp b/packages/thyra/core/src/interfaces/operator_solve/fundamental/Thyra_LinearOpWithSolveBase_decl.hpp index ffd7e8b542e9..6049e827bb02 100644 --- a/packages/thyra/core/src/interfaces/operator_solve/fundamental/Thyra_LinearOpWithSolveBase_decl.hpp +++ b/packages/thyra/core/src/interfaces/operator_solve/fundamental/Thyra_LinearOpWithSolveBase_decl.hpp @@ -279,6 +279,8 @@ class LinearOpWithSolveBase { public: + using scalar_type = Scalar; + /** \name Public interface funtions. */ //@{