Trilinos: use Tpetra without enabling Epetra

[mayrmt]

Current status

When configuring with DEAL_II_WITH_TRILINOS, it is automatically assumed, that Trilinos enables the Epetra linear algebra stack. In addition, one can enable the Tpetra stack via DEAL_II_TRILINOS_WITH_TPETRA (optional).

Epetra is always enabled. Hence, one cannot use Tetra only.

Problem description

In applications, one can only use one linear algebra stack. When moving to Tpetra (the more modern stack), why should one be forced to provide Trilinos with both Epetra and Tpetra?

Possible solution

At first glance, I see two options:

1. Force user to select one stack in the cmake command. CMake logic can then check, if one asks for with exactly one linear algebra stack. Error out otherwise.
2. Move to Xpetra, the Trilinos-internal wrapper to switch between Epetra and Tpetra.

As Xpetra is just a wrapper, I'd recommend the first option, though I wanted to mention the second option for completeness.

I can't estimate the required changes to realize this. I have done the CMake part for option 1 locally (and can push if desired), but haven't adapted the source code yet. There are quite many places where one would need to distinguish between Epetra and Tpetra. I'm wondering if there's a cleaner and more clever way than just scattering #ifdef DEAL_II_TRILINOS_WITH_EPETRA / #ifdef DEAL_II_TRILINOS_WITH_TPETRA all over the code base.

[Rombur]

@mayrmt Yes, that something we should change. Epetra is not optional for historical reason. Another problem is that right now, we don't have wrappers for the entire Tpetra stack (Belos, Amesos2, etc.).
As for the solution, I think option 1 is the way to go. My understanding is that if the Trilinos developers could kill Xpetra, they would.

We will be happy to take any patch that get us closer to have full support for Tpetra.

[masterleinad]

@mayrmt How much success have you had using Tpetra with deal.II? IIRC we don't have support for Tpetra matrices yet.

[bangerth]

I think that's the key question: all of the TrilinosWrappers classes build on Epetra, but we don't have a complete replacement for these classes that would build on Tpetra. We would like to have such replacements (in which case we would just remove the Epetra interfaces) but we lack the manpower to write all of these interfaces -- help is of course always welcome!

[mayrmt]

@Rombur Yes, if Trilinos could get rid of Xpetra, they would. So, I agree that the user needs to activate either the Epetra or Tpetra stack (with the other one being disabled).

@masterleinad @bangerth So far, I have adapted CMake to force the user to either choose Epetra or Tpetra. I'm trying to figure out a path forward to tackle Tpetra-related changes.

I'm now guarding Epetra code with #ifdef DEAL_II_TRILINOS_WITH_EPETRA, such that it actually could be disabled. Progress is slow, but hopefully steady.

At the moment, I'm struggling with Tpetra template parameters: I need to define for example Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> (using template parameter names from Trilinos). What are deal.II's guidelines on choosing them?

[bangerth]

@mayrmt Thanks for working on this!

As for template parameters: If you just want to convert the existing Epetra wrappers, then think Scalar=double, LocalOrdinal=unsigned int and GlobalOrdinal=types::global_dof_index are appropriate, as would be setting both to types::global_dof_index. Whether this choice is acceptable I don't know: types::global_dof_index is an unsigned data type. If that doesn't work for Tpetra, then choosing the corresponding signed integer data types would work as well, i.e., LocalOrdinal=int, and GlobalOrdinal=std::make_signed<types::global_dof_index>::type.

For Node, is there a default that corresponds to just CPU without any additional parallelization or offloading?

[Rombur]

There are plans to get ride of the GlobalOrdinal templates parameters trilinos/Trilinos#2548 This will probably take a long time before it gets done but it might be a good idea to already pick GlobalOrdinal=int64_t

[masterleinad]

For Node, is there a default that corresponds to just CPU without any additional parallelization or offloading?

Kokkos::Compat::KokkosSerialWrapperNode

[jpthiele]

I think that's the key question: all of the TrilinosWrappers classes build on Epetra, but we don't have a complete replacement for these classes that would build on Tpetra. We would like to have such replacements (in which case we would just remove the Epetra interfaces) but we lack the manpower to write all of these interfaces -- help is of course always welcome!

I'd also be interested as I would like to use ShyLU_DD's preconditioners for Krylow methods.
But I can't get the Epetra->Xpetra conversion to work in Trilinos 13.2-
Is there an overview/checklist of what would have to be done?
For just the basic LA (Block) SparseMatrix and Vector come to mind with matching methods for AffineConstraints.

Additionally I saw that Tpetra has NUMA optimizations which would be nice for HPC