Towards Tpetra support

[gassmoeller]

With the changes in #6985 all ASPECT tests now finish the tester system that we also use for the deal.II development version when using the Tpetra vectors and solvers. I have also checked that the solkz and solcx results are similar in accuracy, which they are and the overall test output differences between Epetra and Tpetra solutions are reasonably small (see https://github.com/gassmoeller/aspect/tree/tpetra_test_differences). However, solver performance is another matter, sometimes the Tpetra Stokes solver needs many more iterations than the Epetra based solver (in particular for the Stokes equation, advection is almost equal), I have not checked runtime or scalability yet.

But by now it could be safe to enable using Tpetra, if the user requests it, or the user builds deal.II without Epetra (essentially equal to using Trilinos 17). But we could also keep this PR open for a while. I think eventually someone will have to spend a bit of time with the Stokes solver parameters and check if we should adjust any for Tpetra (similar to the experiments in #234).

[bangerth]

The OS X tester says:

-- Found deal.II version 9.6.2 at '/Users/jenkins/candi-9.6.2-r1/deal.II-v9.6.2/lib/cmake/deal.II'
CMake Error at CMakeLists.txt:131 (message):


  -- deal.II was built without support for Trilinos Epetra and without
  support for Trilinos Tpetra.  This is not supported.



-- Using ASPECT_USE_TPETRA = 'ON'
CMake Error at CMakeLists.txt:145 (message):


  -- deal.II was built without support for Tpetra!



CMake Error at CMakeLists.txt:152 (message):


  -- Tpetra support is only available for deal.II versions greater or equal
  to 9.8.0!

We may simply not have exported DEAL_II_TRILINOS_WITH_EPETRA in 9.6 yet. I think you need to catch this somewhere.

[bangerth]

Other than that, if I read your diffs correctly, you have changes such as

-   Solving Stokes system (AMG)... 53+0 iterations.
+   Solving Stokes system (AMG)... 799+0 iterations.

So this has gotten vastly worse. We use our own dealii::SolverGMRES for this, so that only leaves two possibilities: (i) the MueLu AMG is vastly worse than the ML AMG; (ii) there are bugs in matrix or vector operations. This is unrelated to the current patch, but it's probably still worth thinking about which it may be.

[gassmoeller]

I fixed the support for older deal.II versions, the failing tests are unrelated.

And I agree there are quite a number of tests where in particular the Stokes solver seems to need many more iterations than with Epetra and I was thinking in the same directions. The most obvious candidate to me is that we no longer remove the constant modes from the AMG preconditioner, but a bug is of course also possible. I dont know the reasoning behind the removed option for constant modes in MueLu, maybe it was thought to be no longer necessary. I will need to take a closer look at the MueLu documentation, there is some mention of the near nullspace, but it is not explained in much detail. Maybe we can emulate similar behavior to ML with some additional settings.

[gassmoeller]

I confirmed in at least one test that the support for constant modes (or nullspace removal) makes a huge difference for the solver performance of the MueLu preconditioner. This may not be the whole story (I have seen melt tests without periodicity suffer from worse iteration numbers as well), but a part of the problem.

The test crustal_formation_particles.prm (90 degree spherical shell with periodic boundaries in lateral direction) has the following output with Epetra:

Number of active cells: 768 (on 5 levels)
Number of degrees of freedom: 16,838 (6,402+833+3,201+3,201+3,201)

*** Timestep 0:  t=0 years, dt=0 years
   Solving temperature system... 0 iterations.
   Rebuilding Stokes preconditioner...
   Solving Stokes system (AMG)... 37+0 iterations.

Number of active cells: 1,488 (on 6 levels)
Number of degrees of freedom: 33,549 (12,758+1,654+6,379+6,379+6,379)

*** Timestep 0:  t=0 years, dt=0 years
   Solving temperature system... 0 iterations.
   Advecting particles...  done.
   Rebuilding Stokes preconditioner...
   Solving Stokes system (AMG)... 40+0 iterations.

   Postprocessing:
     Writing graphical output:  output-crustal_formation_particles/solution/solution-00000
     Compositions min/max/mass: 0/0.25/1.976e+11 // 0/0.5412/3.053e+11
     Writing particle output:   output-crustal_formation_particles/particles/particles-00000

*** Timestep 1:  t=1e+06 years, dt=1e+06 years
   Solving temperature system... 8 iterations.
   Advecting particles...  done.
   Solving Stokes system (AMG)... 29+0 iterations.

   Postprocessing:
     Writing graphical output:  output-crustal_formation_particles/solution/solution-00001
     Compositions min/max/mass: 0/0.814/2.217e+11 // 0/0.9645/1.127e+12
     Writing particle output:   output-crustal_formation_particles/particles/particles-00001

Termination requested by criterion: end time

With the current version of MueLu/Tpetra:

Number of active cells: 768 (on 5 levels)
Number of degrees of freedom: 16,838 (6,402+833+3,201+3,201+3,201)

*** Timestep 0:  t=0 years, dt=0 years
   Solving temperature system... 0 iterations.
   Rebuilding Stokes preconditioner...
   Solving Stokes system (AMG)... 73+0 iterations.

Number of active cells: 1,488 (on 6 levels)
Number of degrees of freedom: 33,549 (12,758+1,654+6,379+6,379+6,379)

*** Timestep 0:  t=0 years, dt=0 years
   Solving temperature system... 0 iterations.
   Advecting particles...  done.
   Rebuilding Stokes preconditioner...
   Solving Stokes system (AMG)... 118+0 iterations.

   Postprocessing:
     Writing graphical output:  output/solution/solution-00000
     Compositions min/max/mass: 0/0.25/1.976e+11 // 0/0.5412/3.053e+11
     Writing particle output:   output/particles/particles-00000

*** Timestep 1:  t=1e+06 years, dt=1e+06 years
   Solving temperature system... 8 iterations.
   Advecting particles...  done.
   Solving Stokes system (AMG)... 75+0 iterations.

   Postprocessing:
     Writing graphical output:  output/solution/solution-00001
     Compositions min/max/mass: 0/0.814/2.217e+11 // 0/0.9645/1.127e+12
     Writing particle output:   output/particles/particles-00001

Termination requested by criterion: end time

With MueLu/Tpetra if I remove the periodic boundaries and replace them with tangential boundaries (flow field changes a bit, but not dramatically):

Number of active cells: 768 (on 5 levels)
Number of degrees of freedom: 16,838 (6,402+833+3,201+3,201+3,201)

*** Timestep 0:  t=0 years, dt=0 years
   Solving temperature system... 0 iterations.
   Solving Stokes system (GMG)... 12+0 iterations.

Number of active cells: 1,488 (on 6 levels)
Number of degrees of freedom: 33,549 (12,758+1,654+6,379+6,379+6,379)

*** Timestep 0:  t=0 years, dt=0 years
   Solving temperature system... 0 iterations.
   Advecting particles...  done.
   Solving Stokes system (GMG)... 12+0 iterations.

   Postprocessing:
     Writing graphical output:  output/solution/solution-00000
     Compositions min/max/mass: 0/0.25/1.976e+11 // 0/0.5412/3.053e+11
     Writing particle output:   output/particles/particles-00000

*** Timestep 1:  t=1e+06 years, dt=1e+06 years
   Solving temperature system... 8 iterations.
   Advecting particles...  done.
   Solving Stokes system (GMG)... 12+0 iterations.

   Postprocessing:
     Writing graphical output:  output/solution/solution-00001
     Compositions min/max/mass: 0/0.8781/2.292e+11 // 0/0.9645/1.131e+12
     Writing particle output:   output/particles/particles-00001

Termination requested by criterion: end time

I will need to look closer at how to implement this properly in deal.II. It will probably require porting the function PreconditionAMG::AdditionalData::set_operator_null_space to PreconditionAMGMueLu<Number, MemorySpace>::initialize and assigning the computed nullspace to the MueLu parameters, maybe as follows:

      coarseParams.set("fix nullspace", true);
      coarseParams.set("Nullspace", nullspace_vector);
      parameter_list.set("nullspace: calculate rotations", true);

This will likely be a bigger undertaking. Do you think we need it for this PR? I had mentioned in the changelog entry that performance may not be as good as the ML preconditioner.

[gassmoeller]

With the changes in dealii/dealii#28453 I can now activate the extraction of constant modes for ASPECT in this PR as well and it significantly improves the performance compared to the results I posted above.

I have run the full testsuite with Tpetra MueLu in this commit.

For example the test I compared above crustal_formation_with_particles.prm:

Epetra/ML:

   Solving Stokes system (AMG)... 37+0 iterations.

Tpetra/MueLu without nullspace removal:

   Solving Stokes system (AMG)... 73+0 iterations.

Tpetra/MueLu with nullspace removal:

   Solving Stokes system (AMG)... 52+0 iterations.

This one is still worse than before, but for most tests iteration count is actually within +/- 10% (and quite a few tests, particularly the ones with two-phase flow actually perform better than before).

The one you pointed out Wolfgang (adiabatic_heating_with_melt) is back to normal:
Epetra/ML: Solving Stokes system (AMG)... 53+0 iterations.
Tpetra/ML/no constant modes extracted: Solving Stokes system (AMG)... 799+0 iterations.
Tpetra/ML/constant modes extracted: Solving Stokes system (AMG)... 59+0 iterations.

These tests were not run in our docker tester container so there can be other changes as well, but I think this is close enough in performance that we can say we support Tpetra/MueLu in ASPECT (with deal.II 9.8 and when dealii/dealii#28453 is merged).