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 <pybind11/stl.h>
 +include <pybind11/stl_bind.h>
+
+# 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 <Teuchos_RCP.hpp>
-+custom_shared Teuchos::RCP
+
 +include_for_namespace Teuchos <PyTrilinos2_Teuchos_Custom.hpp>
 +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<void>
+-class Teuchos::ArrayRCP<const void>
+-class Teuchos::Array<bool>
+-class Teuchos::Ptr<int>
+-class Teuchos::Ptr<long>
+-class Teuchos::Ptr<long int>
+-class Teuchos::Ptr<long long>
+-class Teuchos::Ptr<unsigned long>
+-class Teuchos::Ptr<long long int>
+-class Teuchos::Ptr<long unsigned int>
+-class Teuchos::Ptr<bool>
+-class Teuchos::Ptr<double>
+
++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 <PyTrilinos2_Tpetra_Custom.hpp>
-+add_on_binder Tpetra::CrsGraph define_CrsGraph_member_functions
-+include <PyTrilinos2_Tpetra_Types.hpp>
-+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 <Cuda/Kokkos_Cuda_Instance.hpp>
@@ -91,88 +145,71 @@
 -include <OpenMP/Kokkos_OpenMP_Instance.hpp>
 -include <impl/Kokkos_HostThreadTeam.hpp>
 -include <impl/Kokkos_InitializationSettings.hpp>
--class Teuchos::ArrayView<ompi_status_public_t>
--class Teuchos::ArrayView<const ompi_status_public_t>
--class Teuchos::Serializer
--class Teuchos::PromotionTraits
--class Teuchos::CommandLineProcessor
--class Teuchos::Ptr<Teuchos::RCP<Teuchos::CommRequest<int> >>
--class Teuchos::ArrayView<Teuchos::RCP<Teuchos::CommRequest<int> >>
--class Teuchos::TimeMonitor
+-include <View/Kokkos_ViewLegacy.hpp>
+-include <View/Kokkos_ViewTraits.hpp>
+-include <View/Kokkos_ViewMapping.hpp>
+-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 <PyTrilinos2_Tpetra_Custom.hpp>
++include <PyTrilinos2_Tpetra_Types.hpp>
++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<int>
--function Teuchos::getRawMpiComm
--class Teuchos::OpaqueWrapper<ompi_communicator_t *>
--class Teuchos::OpaqueWrapper<ompi_errhandler_t *>
--function Teuchos::Details::setMpiReductionOp
--function Teuchos::Details::getMpiOpForEReductionType
--function Tpetra::Details::extractMpiCommFromTeuchos
--class Teuchos::OpaqueWrapper<ompi_errhandler_t *>
--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 <Stratimikos_MueLuHelpers.hpp>
+
+###############################################################################
+# Thyra
++include <Thyra_VectorSpaceFactoryBase.hpp>
+-class Thyra::ProductVectorSpace
+-class Thyra::RowStatLinearOpBase
+-class Thyra::ScaledLinearOpBase
+-class Thyra::LinearOpSourceBase
++include_for_namespace Thyra <PyTrilinos2_Thyra_Custom.hpp>
++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 <PyTrilinos2_config.hpp>
 #include <PyTrilinos2_Teuchos_ETI.hpp>
 #include <PyTrilinos2_Tpetra_ETI.hpp>
+#include <PyTrilinos2_Thyra_ETI.hpp>
+#include <PyTrilinos2_Stratimikos_ETI.hpp>
+#ifdef HAVE_PYTRILINOS2_MUELU
 #include <PyTrilinos2_MueLu_ETI.hpp>
+#endif
+#include <PyTrilinos2_stream.hpp>
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 <PyTrilinos2_config.hpp>
+#include <Stratimikos_LinearSolverBuilder.hpp>
+#ifdef HAVE_PYTRILINOS2_MUELU
+#include <Stratimikos_MueLuHelpers.hpp>
+#endif
+
+
+#define BINDER_STRATIMIKOS_LINEARSOLVERBUILDER_INSTANT(SC) \
+  template class Stratimikos::LinearSolverBuilder<SC>;
+
+#define BINDER_STRATIMIKOS_MUELU_INSTANT(SCALAR,LO,GO,NO) \
+  template void enableMueLu<SCALAR, LO, GO, NO>(Stratimikos::LinearSolverBuilder<SCALAR>& builder, const std::string& stratName); \
+  template void enableMueLuRefMaxwell<SCALAR, LO, GO, NO>(Stratimikos::LinearSolverBuilder<SCALAR>& builder, const std::string& stratName); \
+  template void enableMueLuMaxwell1<SCALAR, LO, GO, NO>(Stratimikos::LinearSolverBuilder<SCALAR>& builder, const std::string& stratName);
+
+
+namespace Stratimikos {
+
+  template <typename T>
+    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 <Thyra_LinearOpWithSolveBase_decl.hpp>
+
+template<typename T>
+void define_solve(T cl) {
+  using SCALAR = typename T::type::scalar_type;
+
+  cl.def("solve",[](Teuchos::RCP<Thyra::LinearOpWithSolveBase<SCALAR> > &m,
+                    const Thyra::EOpTransp A_trans,
+                    const Thyra::MultiVectorBase<SCALAR> &B,
+                    const Teuchos::RCP<Thyra::MultiVectorBase<SCALAR> > &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 <Thyra_TpetraThyraWrappers.hpp>
+#include <Thyra_TpetraEuclideanScalarProd.hpp>
+#include <Thyra_TpetraVectorSpace.hpp>
+#include <Thyra_TpetraMultiVector.hpp>
+#include <Thyra_TpetraVector.hpp>
+#include <Thyra_TpetraLinearOp.hpp>
+
+#include <Thyra_LinearOpWithSolveBase.hpp>
+#include <Thyra_LinearSolverBuilderBase.hpp>
+#include <Thyra_VectorBase.hpp>
+#include <Thyra_LinearOpWithSolveFactoryHelpers.hpp>
+#include <Thyra_SolveSupportTypes.hpp>
+#include <Thyra_OperatorVectorTypes.hpp>
+
+
+#define BINDER_THYRA_WRAPPERS_INSTANT(SC,LO,GO,NO) \
+  template Teuchos::RCP<Thyra::VectorBase<SC> > createVector<SC, LO, GO, NO>(const Teuchos::RCP<Tpetra::Vector<SC, LO, GO, NO> > &tpetraVector, \
+                                                                             const Teuchos::RCP<const Thyra::VectorSpaceBase<SC> > space); \
+  template Teuchos::RCP<Thyra::MultiVectorBase<SC> > createMultiVector<SC, LO, GO, NO>(const Teuchos::RCP<Tpetra::MultiVector<SC, LO, GO, NO> > &tpetraMultiVector, \
+                                                                                       const Teuchos::RCP<const Thyra::VectorSpaceBase<SC> > rangeSpace, \
+                                                                                       const Teuchos::RCP<const Thyra::VectorSpaceBase<SC> > domainSpace);
+
+#define BINDER_THYRA_EUCLIDEANSCALARPROD_INSTANT(SC,LO,GO,NO) \
+  template class Thyra::TpetraEuclideanScalarProd<SC, LO, GO, NO>;
+
+#define BINDER_THYRA_VECTORSPACE_INSTANT(SC,LO,GO,NO) \
+  template class Thyra::TpetraVectorSpace<SC, LO, GO, NO>;
+
+#define BINDER_THYRA_MULTIVECTOR_INSTANT(SC,LO,GO,NO) \
+  template class Thyra::TpetraMultiVector<SC, LO, GO, NO>;
+
+#define BINDER_THYRA_VECTOR_INSTANT(SC, LO, GO, NO) \
+  template class Thyra::TpetraVector<SC, LO, GO, NO>;
+
+#define BINDER_THYRA_LINEAROP_INSTANT(SC,LO,GO,NO) \
+  template class Thyra::TpetraLinearOp<SC, LO, GO, NO>; \
+  template Teuchos::RCP<Thyra::TpetraLinearOp<SC, LO, GO, NO> > tpetraLinearOp(const Teuchos::RCP<const Thyra::VectorSpaceBase<SC> > &rangeSpace, \
+                                                                               const Teuchos::RCP<const Thyra::VectorSpaceBase<SC> > &domainSpace, \
+                                                                               const Teuchos::RCP<Tpetra::Operator<SC, LO, GO, NO> > &tpetraOperator); \
+  template Teuchos::RCP<const Thyra::TpetraLinearOp<SC, LO, GO, NO> > constTpetraLinearOp(const Teuchos::RCP<const Thyra::VectorSpaceBase<SC> > &rangeSpace, \
+                                                                                          const Teuchos::RCP<const Thyra::VectorSpaceBase<SC> > &domainSpace, \
+                                                                                          const Teuchos::RCP<const Tpetra::Operator<SC, LO, GO, NO> > &tpetraOperator);
+
+#define BINDER_THYRA_LINEARSOLVER_INSTANT(SC) \
+  template class Thyra::LinearOpWithSolveFactoryBase<SC>; \
+  template Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<SC> > Thyra::createLinearSolveStrategy(const Thyra::LinearSolverBuilderBase<SC> &linearSolverBuilder, const std::string &linearSolveStrategyName); \
+  template Teuchos::RCP<Thyra::LinearOpWithSolveBase<SC> > linearOpWithSolve(const Thyra::LinearOpWithSolveFactoryBase<SC> &lowsFactory, \
+                                                                             const Teuchos::RCP<const Thyra::LinearOpBase<SC> > &fwdOp, \
+                                                                             const Thyra::ESupportSolveUse supportSolveUse); \
+  template struct Thyra::SolveCriteria<SC>; \
+  template struct Thyra::SolveStatus<SC>;
+
+#define BINDER_TPETRA_OPERATOR_INSTANT(SC,LO,GO,NO) \
+   template class Tpetra::Operator<SC, LO, GO, NO>;
+
+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 <typename T>
+    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 <typename T>
 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 <iostream>
+#include <fstream>
+#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<std::ostream> getCout() {
+    Teuchos::RCP<std::ostream> 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/muelu/adapters/CMakeLists.txt b/packages/muelu/adapters/CMakeLists.txt
index f6498df28a69..7a4314b74552 100644
--- a/packages/muelu/adapters/CMakeLists.txt
+++ b/packages/muelu/adapters/CMakeLists.txt
@@ -112,6 +112,7 @@ IF (${PACKAGE_NAME}_ENABLE_Stratimikos)
   IF (${PACKAGE_NAME}_ENABLE_Stratimikos AND ${PACKAGE_NAME}_ENABLE_ThyraTpetraAdapters)
     TRIBITS_SET_AND_INC_DIRS(DIR ${CMAKE_CURRENT_SOURCE_DIR}/stratimikos)
     IF (${PACKAGE_NAME}_ENABLE_Teko AND ${PACKAGE_NAME}_ENABLE_Experimental)
+      TRIBITS_INCLUDE_DIRECTORIES(${DIR}/../../research/q2q1)
       APPEND_SET(HEADERS
         ${DIR}/Thyra_MueLuTpetraQ2Q1PreconditionerFactory_decl.hpp
         ${DIR}/Thyra_MueLuTpetraQ2Q1PreconditionerFactory_def.hpp
diff --git a/packages/muelu/adapters/stratimikos/Thyra_MueLuTpetraQ2Q1PreconditionerFactory_def.hpp b/packages/muelu/adapters/stratimikos/Thyra_MueLuTpetraQ2Q1PreconditionerFactory_def.hpp
index 46d9ac613009..7892138358fb 100644
--- a/packages/muelu/adapters/stratimikos/Thyra_MueLuTpetraQ2Q1PreconditionerFactory_def.hpp
+++ b/packages/muelu/adapters/stratimikos/Thyra_MueLuTpetraQ2Q1PreconditionerFactory_def.hpp
@@ -36,8 +36,8 @@
 
 #include "MueLu.hpp"
 
-#include "../../research/q2q1/MueLu_Q2Q1PFactory.hpp"
-#include "../../research/q2q1/MueLu_Q2Q1uPFactory.hpp"
+#include "MueLu_Q2Q1PFactory.hpp"
+#include "MueLu_Q2Q1uPFactory.hpp"
 
 #include "MueLu_AmalgamationFactory.hpp"
 #include "MueLu_BaseClass.hpp"
diff --git a/packages/muelu/adapters/tpetra/MueLu_TpetraOperatorAsRowMatrix.hpp b/packages/muelu/adapters/tpetra/MueLu_TpetraOperatorAsRowMatrix.hpp
index 379e830ca577..a096f3973080 100644
--- a/packages/muelu/adapters/tpetra/MueLu_TpetraOperatorAsRowMatrix.hpp
+++ b/packages/muelu/adapters/tpetra/MueLu_TpetraOperatorAsRowMatrix.hpp
@@ -6,8 +6,12 @@
 // SPDX-License-Identifier: BSD-3-Clause
 // *****************************************************************************
 // @HEADER
+#ifndef MUELU_TPETRAOPERATORASROWMATRIX_HPP
+#define MUELU_TPETRAOPERATORASROWMATRIX_HPP
 
+#include <Teuchos_RCP.hpp>
 #include <Tpetra_RowMatrix.hpp>
+#include <MueLu_Exceptions.hpp>
 
 namespace MueLu {
 
@@ -244,3 +248,5 @@ class TpetraOperatorAsRowMatrix : public Tpetra::RowMatrix<Scalar, LocalOrdinal,
 };
 
 }  // namespace MueLu
+
+#endif
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. */
   //@{