Porting tutorials to PETSc from Trilinos

[Franco Milicchio]

Dear all,

I am trying to port the Trilinos examples, for instance Step-31, to PETSc. I am finding some difficulties, though, in the block matrix class, it seems a little bit different from Trilinos. In this function:

  template <int dim>

  void BoussinesqFlowProblem<dim>::solve()

  {

    std::cout << "   Solving..." << std::endl;

    {

        const LinearSolvers::InverseMatrix<PETScWrappers::SparseMatrix,

                                         PETScWrappers::PreconditionICC>

        mp_inverse(stokes_preconditioner_matrix.block(1, 1),

                   *Mp_preconditioner);

the LinearSolvers::InverseMatrix object can't be initialized

error: no matching constructor for initialization of 'const LinearSolvers::InverseMatrix<PETScWrappers::SparseMatrix, PETScWrappers::PreconditionICC>'

note: candidate constructor not viable: no known conversion from 'dealii::BlockMatrixBase<dealii::PETScWrappers::MPI::SparseMatrix>::BlockType' (aka 'dealii::PETScWrappers::MPI::SparseMatrix') to 'const dealii::SparseMatrix<double>' for 1st argument

As I didn't find a 1:1 correspondence between the two namespaces, I've resorted to remove Trilinos classes as best as I could, for instance, TrilinosWrappers::BlockSparseMatrix exists but in PETSc there is only PETScWrappers::MPI::BlockSparseMatrix that I know of. And the two seem to be not exactly equivalent.

The complete Trilinos/PETSc variables are the following:

// ORIGINAL TRILINOS

TrilinosWrappers::BlockSparseMatrix stokes_matrix;

TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix;

TrilinosWrappers::MPI::BlockVector stokes_solution;

TrilinosWrappers::MPI::BlockVector old_stokes_solution;

TrilinosWrappers::MPI::BlockVector stokes_rhs;

TrilinosWrappers::SparseMatrix temperature_mass_matrix;

TrilinosWrappers::SparseMatrix temperature_stiffness_matrix;

TrilinosWrappers::SparseMatrix temperature_matrix;

TrilinosWrappers::MPI::Vector temperature_solution;

TrilinosWrappers::MPI::Vector old_temperature_solution;

TrilinosWrappers::MPI::Vector old_old_temperature_solution;

TrilinosWrappers::MPI::Vector temperature_rhs;

std::shared_ptr<TrilinosWrappers::PreconditionAMG> Amg_preconditioner;

std::shared_ptr<TrilinosWrappers::PreconditionIC>  Mp_preconditioner;

// NEW VARIABLES WITH PETSC

PETScWrappers::MPI::BlockSparseMatrix stokes_matrix;

PETScWrappers::MPI::BlockSparseMatrix stokes_preconditioner_matrix;

PETScWrappers::MPI::BlockVector stokes_solution;

PETScWrappers::MPI::BlockVector old_stokes_solution;

PETScWrappers::MPI::BlockVector stokes_rhs;

PETScWrappers::SparseMatrix temperature_mass_matrix;

PETScWrappers::SparseMatrix temperature_stiffness_matrix;

PETScWrappers::SparseMatrix temperature_matrix;

PETScWrappers::MPI::Vector temperature_solution;

PETScWrappers::MPI::Vector old_temperature_solution;

PETScWrappers::MPI::Vector old_old_temperature_solution;

PETScWrappers::MPI::Vector temperature_rhs;

std::shared_ptr<PETScWrappers::PreconditionILU> Amg_preconditioner;

std::shared_ptr<PETScWrappers::PreconditionICC>  Mp_preconditioner;

What classes am I missing here?

Thanks!

[Daniel Arndt]

Franco,

 

[...]

As I didn't find a 1:1 correspondence between the two namespaces, I've resorted to remove Trilinos classes as best as I could, for instance, TrilinosWrappers::BlockSparseMatrix exists but in PETSc there is only PETScWrappers::MPI::BlockSparseMatrix that I know of. And the two seem to be not exactly equivalent.

step-40 tries to allow switching between Trilinos and PETSc as easy as possible by using the deal.II/lac/generic_linear_algebra.h header and the

defining the LA::MPI namespace.

 

[...]

std::shared_ptr<PETScWrappers::PreconditionILU> Amg_preconditioner;

 

This looks inappropriate. Probably it should rather be PETScWrappers::PreconditionBoomerAMG;

Best,

Daniel

[Franco Milicchio]

step-40 tries to allow switching between Trilinos and PETSc as easy as possible by using the deal.II/lac/generic_linear_algebra.h header and the

defining the LA::MPI namespace.

Thanks for your answer, Daniel. Even with that example I find it a little harder than I expected to use classes from PETSc instead of Trilinos. For instance:

        const LinearSolvers::InverseMatrix<PETScWrappers::SparseMatrix, PETScWrappers::PreconditionICC>

            mp_inverse(stokes_preconditioner_matrix.block(1, 1), *Mp_preconditioner);

Cannot compile, I am trying to use a static_cast<PETScWrappers::SparseMatrix>(stokes_preconditioner_matrix.block(1, 1)) now, but I don't think it will be a good idea.

 

Another weak point in my porting is

    template <class PreconditionerTypeA, class PreconditionerTypeMp>

    BlockSchurPreconditioner<PreconditionerTypeA, PreconditionerTypeMp>::

      BlockSchurPreconditioner(

        const PETScWrappers::MPI::BlockSparseMatrix &S,

        const InverseMatrix<PETScWrappers::SparseMatrix,

                            PreconditionerTypeMp> &Mpinv,

        const PreconditionerTypeA &                Apreconditioner)

      : stokes_matrix(&S)

      , m_inverse(&Mpinv)

      , a_preconditioner(Apreconditioner)

      , tmp(complete_index_set(stokes_matrix->block(1, 1).m()))

    {}

where tmp (a mutable PETScWrappers::MPI::Vector) cannot be initialized. Adding MPI_COMM_WORLD make it compile, though.

Again with matrices, I have made the following changes:

        //temperature_matrix.copy_from(temperature_mass_matrix);

        temperature_matrix = 0;

        temperature_matrix.add(1.0, temperature_mass_matrix);

I hope the change works...

Can I ask why there are some discrepancies between those two wrappers? Apart from an obvious one, that they're different libraries, but uniformity may help a lot IMO.

 

[...]

std::shared_ptr<PETScWrappers::PreconditionILU> Amg_preconditioner;

 

This looks inappropriate. Probably it should rather be PETScWrappers::PreconditionBoomerAMG;

Yes, I know :) I am first trying to make it compile, then I will try to use a good preconditioner (Boomer requires Hypre, and it needs MPI, so right now it is not an option).

Thanks!

    Franco

[Bruno Turcksin]

Franco,

On Tuesday, July 16, 2019 at 3:54:34 AM UTC-4, Franco Milicchio wrote:

Can I ask why there are some discrepancies between those two wrappers? Apart from an obvious one, that they're different libraries, but uniformity may help a lot IMO.

Because the wrappers have been written over time by different people who were not aware of how exactly the other wrapper works (most people use either Trilinos or PETSc not both). We are trying to get the interfaces to be more similar but it is very difficult because any change in the interface will break many codes.

Best,

Bruno

[Franco Milicchio]

Thanks for your answer. By the way, would a wrapper for the Eigen library be appealing to anyone?

And back to my lowly question :)

I was able to port the Step 31 code to plain deal.II without PETSc, but I have an error at runtime. At the moment I removed the preconditioners and just used a PreconditionIdentity to make it compile and link, I will deal with the right preconditioner later.

Apparently I have initialized something wrong, as the assignment to the matrix

    if (rebuild_stokes_matrix == true)

      stokes_matrix = 0;

makes the program die. 

Number of active cells: 256 (on 5 levels)

Number of degrees of freedom: 3556 (n_u 2178 + n_p 289 + n_T 1089)

Timestep 0:  t=0

   Assembling...stokes_matrix 2467 x 2467 blocks 2 x 2

--------------------------------------------------------

An error occurred in line <382> of file </Users/sensei/Downloads/src/dealii-9.1.0/include/deal.II/base/smartpointer.h> in function

    T *dealii::SmartPointer<const dealii::SparsityPattern, dealii::SparseMatrix<double> >::operator->() const [T = const dealii::SparsityPattern, P = dealii::SparseMatrix<double>]

The violated condition was: 

    pointed_to_object_is_alive

Additional information: 

    The object pointed to is not valid anymore.

Stacktrace:

-----------

#0  2   libdeal_II.g.9.1.0.dylib            0x00000001174e003e _ZNK6dealii12SmartPointerIKNS_15SparsityPatternENS_12SparseMatrixIdEEEptEv + 206: 2   libdeal_II.g.9.1.0.dylib            0x00000001174e003e _ZNK6dealii12SmartPointerIKNS_15SparsityPatternENS_12SparseMatrixIdEEEptEv 

#1  3   libdeal_II.g.9.1.0.dylib            0x000000011abc65e5 _ZN6dealii12SparseMatrixIdEaSEd + 69: 3   libdeal_II.g.9.1.0.dylib            0x000000011abc65e5 _ZN6dealii12SparseMatrixIdEaSEd 

#2  4   step-31                             0x000000010ada5f3a _ZN6dealii17BlockSparseMatrixIdEaSEd + 90: 4   step-31                             0x000000010ada5f3a _ZN6dealii17BlockSparseMatrixIdEaSEd 

#3  5   step-31                             0x000000010ada0e77 _ZN6Step3121BoussinesqFlowProblemILi2EE22assemble_stokes_systemEv + 247: 5   step-31                             0x000000010ada0e77 _ZN6Step3121BoussinesqFlowProblemILi2EE22assemble_stokes_systemEv 

#4  6   step-31                             0x000000010ad944b8 _ZN6Step3121BoussinesqFlowProblemILi2EE3runEv + 504: 6   step-31                             0x000000010ad944b8 _ZN6Step3121BoussinesqFlowProblemILi2EE3runEv 

#5  7   step-31                             0x000000010ad93ebd main + 125: 7   step-31                             0x000000010ad93ebd main 

#6  8   libdyld.dylib                       0x00007fff6bbd93d5 start + 1: 8   libdyld.dylib                       0x00007fff6bbd93d5 start 

--------------------------------------------------------

I have checked the sizes as you can see, and they seem to be properly initialized. So my guess is that I am doing something wrong with sparsity patterns, but I don't know how:

      BlockSparsityPattern sp(2, 2);

      sp.copy_from(dsp);

      stokes_matrix.reinit(sp); //stokes_matrix.reinit(dsp);

      // ...

      SparsityPattern sp(n_T, n_T);

      sp.copy_from(dsp);

      temperature_matrix.reinit(sp); //temperature_matrix.reinit(dsp);

      temperature_mass_matrix.reinit(sp); //temperature_mass_matrix.reinit(temperature_matrix);

      temperature_stiffness_matrix.reinit(sp); //temperature_stiffness_matrix.reinit(temperature_matrix);

The code I've modified is available in full here: https://gist.github.com/fmilicchio/2dac430d7c071f98b4f66d79d2101b04

Any hints?

Thank you very much for your help!

    Franco

[Bruno Turcksin]

Franco,

Le jeu. 18 juil. 2019 à 10:27, Franco Milicchio
<franco.m...@gmail.com> a écrit :

> Thanks for your answer. By the way, would a wrapper for the Eigen library be appealing to anyone?

I am sure it would be useful!

> --------------------------------------------------------
> An error occurred in line <382> of file </Users/sensei/Downloads/src/dealii-9.1.0/include/deal.II/base/smartpointer.h> in function
> T *dealii::SmartPointer<const dealii::SparsityPattern, dealii::SparseMatrix<double> >::operator->() const [T = const dealii::SparsityPattern, P = dealii::SparseMatrix<double>]
> The violated condition was:
> pointed_to_object_is_alive
> Additional information:
> The object pointed to is not valid anymore.

It looks like the sparsity pattern does not exist anymore. The matrix
object does not keep a copy of the sparsity pattern, so the sparsity
pattern needs to exist as long as the matrix that uses it is alive.

Best,

Bruno

[Franco Milicchio]

> Thanks for your answer. By the way, would a wrapper for the Eigen library be appealing to anyone? 
I am sure it would be useful! 

Good to know, Bruno. I will start looking into those wrappers and I hope some student of mine will be helping (fingers crossed). 

Is there a small documentation that I can read? For instance, I see there are several namespaces like LinearAlgebra, PETSc/Epetra/Trilinos/CUDAWrappers, LinearAlgebraDealii and so on, that might be restructured without compromising old codes, adding some uniformity to classes. 

 

It looks like the sparsity pattern does not exist anymore. The matrix 
object does not keep a copy of the sparsity pattern, so the sparsity 
pattern needs to exist as long as the matrix that uses it is alive. 

Thanks, that did the trick! The problem now moved to the AffineConstraints (which I didn't touch) complaining about dimensions:

--------------------------------------------------------

An error occurred in line <3781> of file </Users/sensei/Downloads/src/dealii-9.1.0/include/deal.II/lac/affine_constraints.templates.h> in function

    void dealii::AffineConstraints<double>::distribute_local_to_global(const FullMatrix<number> &, const Vector<number> &, const std::vector<size_type> &, MatrixType &, VectorType &, bool, std::integral_constant<bool, true>) const [number = double, MatrixType = dealii::BlockSparseMatrix<double>, VectorType = dealii::BlockVector<double>]

The violated condition was: 

    static_cast<typename ::dealii::internal::argument_type<void( typename std::common_type<decltype(global_matrix.m()), decltype(global_vector.size())>::type)>::type>(global_matrix.m()) == static_cast<typename ::dealii::internal::argument_type<void( typename std::common_type<decltype(global_matrix.m()), decltype(global_vector.size())>::type)>::type>(global_vector.size())

Additional information: 

    Dimension 2467 not equal to 0.

Stacktrace:

-----------

#0  2   libdeal_II.g.9.1.0.dylib            0x00000001164eb3a2 _ZNK6dealii17AffineConstraintsIdE26distribute_local_to_globalINS_17BlockSparseMatrixIdEENS_11BlockVectorIdEEEEvRKNS_10FullMatrixIdEERKNS_6VectorIdEERKNSt3__16vectorIjNSF_9allocatorIjEEEERT_RT0_bNSF_17integral_constantIbLb1EEE + 1570: 2   libdeal_II.g.9.1.0.dylib            0x00000001164eb3a2 _ZNK6dealii17AffineConstraintsIdE26distribute_local_to_globalINS_17BlockSparseMatrixIdEENS_11BlockVectorIdEEEEvRKNS_10FullMatrixIdEERKNS_6VectorIdEERKNSt3__16vectorIjNSF_9allocatorIjEEEERT_RT0_bNSF_17integral_constantIbLb1EEE 

#1  3   step-31                             0x0000000108d5b641 _ZN6Step3121BoussinesqFlowProblemILi2EE22assemble_stokes_systemEv + 1953: 3   step-31                             0x0000000108d5b641 _ZN6Step3121BoussinesqFlowProblemILi2EE22assemble_stokes_systemEv 

#2  4   step-31                             0x0000000108d4e468 _ZN6Step3121BoussinesqFlowProblemILi2EE3runEv + 504: 4   step-31                             0x0000000108d4e468 _ZN6Step3121BoussinesqFlowProblemILi2EE3runEv 

#3  5   step-31                             0x0000000108d4de6d main + 125: 5   step-31                             0x0000000108d4de6d main 

#4  6   libdyld.dylib                       0x00007fff6bbd93d5 start + 1: 6   libdyld.dylib                       0x00007fff6bbd93d5 start 

#5  7   ???                                 0x0000000000000001 0x0 + 1: 7   ???                                 0x0000000000000001 0x0 

--------------------------------------------------------

Signal: SIGABRT (signal SIGABRT)

on this line (initial timestep):

        if (rebuild_stokes_matrix == true)

          stokes_constraints.distribute_local_to_global(local_matrix,

                                                        local_rhs,

                                                        local_dof_indices,

                                                        stokes_matrix,

                                                        stokes_rhs);

The stokes_constraints.distribute_local_to_global apparently doesn't like how I implemented the sparsity patterns to build the matrices. All I did here is to add two sparsity patterns, one for the Stokes matrix, and one for the temperature:

      stokes_matrix_sp.reinit(2, 2); // probably unnecessary

      stokes_matrix_sp.copy_from(dsp);

      stokes_matrix.reinit(stokes_matrix_sp); //stokes_matrix.reinit(dsp);

    BlockSparsityPattern stokes_matrix_sp; // the definition

I am probably doing something terribly stupid here...

Thanks for any help!

    Franco

[Bruno Turcksin]

Franco

Le ven. 19 juil. 2019 à 10:04, Franco Milicchio
<franco.m...@gmail.com> a écrit :

> Is there a small documentation that I can read? For instance, I see there are several namespaces like LinearAlgebra, PETSc/Epetra/Trilinos/CUDAWrappers, LinearAlgebraDealii and so on, that might be restructured without compromising old codes, adding some uniformity to classes.

No unfortunately we don't have such things. However, most of the new
effort on linear algebra is in the LinearAlgebra namespace.

> --------------------------------------------------------
> An error occurred in line <3781> of file </Users/sensei/Downloads/src/dealii-9.1.0/include/deal.II/lac/affine_constraints.templates.h> in function
> void dealii::AffineConstraints<double>::distribute_local_to_global(const FullMatrix<number> &, const Vector<number> &, const std::vector<size_type> &, MatrixType &, VectorType &, bool, std::integral_constant<bool, true>) const [number = double, MatrixType = dealii::BlockSparseMatrix<double>, VectorType = dealii::BlockVector<double>]
> The violated condition was:
> static_cast<typename ::dealii::internal::argument_type<void( typename std::common_type<decltype(global_matrix.m()), decltype(global_vector.size())>::type)>::type>(global_matrix.m()) == static_cast<typename ::dealii::internal::argument_type<void( typename std::common_type<decltype(global_matrix.m()), decltype(global_vector.size())>::type)>::type>(global_vector.size())
> Additional information:

I don't know why this error is so unreadable but if you look at line
3781 in include/deal.II/lac/affine_constraints.templates.h you will
see that the assert checks that the global matrix has as many rows has
the global vector. To me it looks like the problem is that the global
vector as a size of zero.

If you forget about the static_cast in the error message just above
you can actually see what happens:
The violated condition was:
global_matrix.m() == global_vector.size()

Best,

Bruno

[Franco Milicchio]

 

No unfortunately we don't have such things. However, most of the new
effort on linear algebra is in the LinearAlgebra namespace.

Thanks, Bruno, I've started looking into that namespace. There are several options to make it work in a coherent way, I'll ask for comments on these solutions as I figured all classes out (and how to use MPI seamlessly, at first I am thinking about tag dispatching to avoid breaking legacy codes).

 

I don't know why this error is so unreadable but if you look at line
3781 in include/deal.II/lac/affine_constraints.templates.h you will
see that the assert checks that the global matrix has as many rows has
the global vector. To me it looks like the problem is that the global
vector as a size of zero.

If you forget about the static_cast in the error message just above
you can actually see what happens:
The violated condition was:
global_matrix.m() == global_vector.size()

I am lost in this, I have made some stupid error here but I don't see how. I am posting the code here:

     https://gist.github.com/fmilicchio/2dac430d7c071f98b4f66d79d2101b04

As far as I can tell, I've initialized every sparsity pattern, matrix, and vector as they were with Trilinos, but I am surely seeing the code I intended to write, not the code I have actually written.

Do you see any stupid error at first sight?

Thank you!

    Franco

 

[Daniel Arndt]

Franco,

If you forget about the static_cast in the error message just above
you can actually see what happens:
The violated condition was:
global_matrix.m() == global_vector.size()

I am lost in this, I have made some stupid error here but I don't see how. I am posting the code here:

     https://gist.github.com/fmilicchio/2dac430d7c071f98b4f66d79d2101b04

As far as I can tell, I've initialized every sparsity pattern, matrix, and vector as they were with Trilinos, but I am surely seeing the code I intended to write, not the code I have actually written.

You are not initializing the stokes vectors appropriately. Running in a debugger should have told you that there was a problem with the size of stokes_rhs. In fact,

stokes_solution.reinit(n_u+n_p)

creates a vector with n_u+n_p empty blocks. You want to do something like

std::vector<types::global_dof_index> block_sizes(2);
block_sizes[0] = n_u;
block_sizes[1] = n_p;
stokes_solution.reinit(block_sizes);

instead. Then, you also need to have a look how to initialize the sparsity pattern for the stokes_preconditioner correctly (if you want to use that matrix at all).

Best,

Daniel

[Franco Milicchio]

Thanks Daniel, now it runs.

Of course it won't converge, lacking preconditiones, but this is for another question.

Thank you!

    Franco