Merge pull request #16625 from jpthiele/tpetra_test_dealii_solvers

Add tests for deal.II solvers with TpetraWrappers classes

include/deal.II/lac/trilinos_tpetra_vector.h

@@ -366,6 +366,24 @@ namespace LinearAlgebra
       reinit(const Vector<Number, MemorySpace> &V,
              const bool                         omit_zeroing_entries = false);
 
+      /**
+       * Swap the contents of this vector and the other vector @p v. One could do
+       * this operation with a temporary variable and copying over the data
+       * elements, but this function is significantly more efficient since it
+       * only swaps the pointers to the data of the two vectors and therefore
+       * does not need to allocate temporary storage and move data around.
+       *
+       * This function is analogous to the @p swap function of all C++
+       * standard containers. Also, there is a global function
+       * <tt>swap(u,v)</tt> that simply calls <tt>u.swap(v)</tt>, again in
+       * analogy to standard functions.
+       *
+       * This function is virtual in order to allow for derived classes to
+       * handle memory separately.
+       */
+      virtual void
+      swap(Vector &v);
+
       /**
        * Extract a range of elements all at once.
        */
@@ -639,6 +657,16 @@ namespace LinearAlgebra
       real_type
       linfty_norm() const;
 
+      /**
+       * Return the square of the l<sub>2</sub>-norm.
+       */
+      real_type
+      norm_sqr() const;
+
+      /**
+       *
+       *
+       */
       /**
        * Performs a combined operation of a vector addition and a subsequent
        * inner product, returning the value of the inner product. In other
@@ -923,6 +951,14 @@ namespace LinearAlgebra
 
     /* ------------------------- Inline functions ---------------------- */
 
+    template <typename Number, typename MemorySpace>
+    inline void
+    swap(Vector<Number, MemorySpace> &u, Vector<Number, MemorySpace> &v)
+    {
+      u.swap(v);
+    }
+
+
     template <typename Number, typename MemorySpace>
     inline bool
     Vector<Number, MemorySpace>::has_ghost_elements() const
@@ -939,6 +975,12 @@ namespace LinearAlgebra
       return compressed;
     }
 
+    template <typename Number, typename MemorySpace>
+    inline void
+    Vector<Number, MemorySpace>::swap(Vector<Number, MemorySpace> &v)
+    {
+      vector.swap(v.vector);
+    }
 
 
     template <typename Number, typename MemorySpace>

include/deal.II/lac/trilinos_tpetra_vector.templates.h

@@ -876,6 +876,15 @@ namespace LinearAlgebra
 
 
 
+    template <typename Number, typename MemorySpace>
+    typename Vector<Number, MemorySpace>::real_type
+    Vector<Number, MemorySpace>::norm_sqr() const
+    {
+      Vector<Number, MemorySpace>::real_type d = l2_norm();
+      return d * d;
+    }
+
+
     template <typename Number, typename MemorySpace>
     Number
     Vector<Number, MemorySpace>::add_and_dot(

tests/trilinos_tpetra/deal_solver_01.cc

@@ -0,0 +1,87 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2004 - 2024 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// test the CG solver using the Trilinos matrix and vector classes
+
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_bicgstab.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/solver_gmres.h>
+#include <deal.II/lac/solver_qmrs.h>
+#include <deal.II/lac/solver_richardson.h>
+#include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
+#include <deal.II/lac/vector_memory.h>
+
+#include <iostream>
+#include <typeinfo>
+
+#include "../tests.h"
+
+#include "../testmatrix.h"
+
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+  deallog << std::setprecision(4);
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(
+    argc, argv, testing_max_num_threads());
+
+
+  {
+    SolverControl control(200, 1.e-3);
+
+    const unsigned int size = 32;
+    unsigned int       dim  = (size - 1) * (size - 1);
+
+    deallog << "Size " << size << " Unknowns " << dim << std::endl;
+
+    // Make matrix
+    FDMatrix               testproblem(size, size);
+    DynamicSparsityPattern csp(dim, dim);
+    testproblem.five_point_structure(csp);
+    LinearAlgebra::TpetraWrappers::SparseMatrix<double> A;
+    A.reinit(csp);
+    testproblem.five_point(A);
+
+    LinearAlgebra::TpetraWrappers::Vector<double> f;
+    f.reinit(complete_index_set(dim), MPI_COMM_WORLD);
+    LinearAlgebra::TpetraWrappers::Vector<double> u;
+    u.reinit(complete_index_set(dim), MPI_COMM_WORLD);
+    f.trilinos_rcp()->putScalar(1.0);
+    A.compress(VectorOperation::insert);
+    f.compress(VectorOperation::insert);
+    u.compress(VectorOperation::insert);
+
+    GrowingVectorMemory<LinearAlgebra::TpetraWrappers::Vector<double>> mem;
+    SolverCG<LinearAlgebra::TpetraWrappers::Vector<double>> solver(control,
+                                                                   mem);
+    PreconditionIdentity                                    preconditioner;
+
+    deallog << "Solver type: " << typeid(solver).name() << std::endl;
+
+    check_solver_within_range(solver.solve(A, u, f, preconditioner),
+                              control.last_step(),
+                              42,
+                              44);
+  }
+}

tests/trilinos_tpetra/deal_solver_01.output

@@ -0,0 +1,4 @@
+
+DEAL::Size 32 Unknowns 961
+DEAL::Solver type: N6dealii8SolverCGINS_13LinearAlgebra14TpetraWrappers6VectorIdNS_11MemorySpace4HostEEEEE
+DEAL::Solver stopped within 42 - 44 iterations

tests/trilinos_tpetra/deal_solver_02.cc

@@ -0,0 +1,88 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2004 - 2024 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// test the BiCGStab solver using the Trilinos matrix and vector classes
+
+
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_bicgstab.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/solver_gmres.h>
+#include <deal.II/lac/solver_qmrs.h>
+#include <deal.II/lac/solver_richardson.h>
+#include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
+#include <deal.II/lac/vector_memory.h>
+
+#include <iostream>
+#include <typeinfo>
+
+#include "../tests.h"
+
+#include "../testmatrix.h"
+
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+  deallog << std::setprecision(4);
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(
+    argc, argv, testing_max_num_threads());
+
+
+  {
+    SolverControl control(200, 1.3e-10);
+
+    const unsigned int size = 32;
+    unsigned int       dim  = (size - 1) * (size - 1);
+
+    deallog << "Size " << size << " Unknowns " << dim << std::endl;
+
+    // Make matrix
+    FDMatrix               testproblem(size, size);
+    DynamicSparsityPattern csp(dim, dim);
+    testproblem.five_point_structure(csp);
+    LinearAlgebra::TpetraWrappers::SparseMatrix<double> A;
+    A.reinit(csp);
+    testproblem.five_point(A);
+
+    LinearAlgebra::TpetraWrappers::Vector<double> f;
+    f.reinit(complete_index_set(dim), MPI_COMM_WORLD);
+    LinearAlgebra::TpetraWrappers::Vector<double> u;
+    u.reinit(complete_index_set(dim), MPI_COMM_WORLD);
+    f.trilinos_rcp()->putScalar(1.0);
+    A.compress(VectorOperation::insert);
+    f.compress(VectorOperation::insert);
+    u.compress(VectorOperation::insert);
+
+    GrowingVectorMemory<LinearAlgebra::TpetraWrappers::Vector<double>> mem;
+    SolverBicgstab<LinearAlgebra::TpetraWrappers::Vector<double>>      solver(
+      control, mem);
+    PreconditionIdentity preconditioner;
+
+    deallog << "Solver type: " << typeid(solver).name() << std::endl;
+
+    check_solver_within_range(solver.solve(A, u, f, preconditioner),
+                              control.last_step(),
+                              49,
+                              51);
+  }
+}

tests/trilinos_tpetra/deal_solver_02.output

@@ -0,0 +1,4 @@
+
+DEAL::Size 32 Unknowns 961
+DEAL::Solver type: N6dealii14SolverBicgstabINS_13LinearAlgebra14TpetraWrappers6VectorIdNS_11MemorySpace4HostEEEEE
+DEAL::Solver stopped within 49 - 51 iterations

tests/trilinos_tpetra/deal_solver_03.cc

@@ -0,0 +1,90 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2004 - 2024 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// test the GMRES solver using the Trilinos matrix and vector classes
+
+
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver.h>
+#include <deal.II/lac/solver_bicgstab.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/solver_gmres.h>
+#include <deal.II/lac/solver_qmrs.h>
+#include <deal.II/lac/solver_richardson.h>
+#include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/vector_memory.h>
+
+#include <iostream>
+#include <typeinfo>
+
+#include "../tests.h"
+
+#include "../testmatrix.h"
+
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+  deallog << std::setprecision(2);
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(
+    argc, argv, testing_max_num_threads());
+
+
+  {
+    SolverControl control(2000, 1.e-3);
+
+    const unsigned int size = 32;
+    unsigned int       dim  = (size - 1) * (size - 1);
+
+    deallog << "Size " << size << " Unknowns " << dim << std::endl;
+
+    // Make matrix
+    FDMatrix               testproblem(size, size);
+    DynamicSparsityPattern csp(dim, dim);
+    testproblem.five_point_structure(csp);
+    LinearAlgebra::TpetraWrappers::SparseMatrix<double> A;
+    A.reinit(csp);
+    testproblem.five_point(A);
+
+    LinearAlgebra::TpetraWrappers::Vector<double> f;
+    f.reinit(complete_index_set(dim), MPI_COMM_WORLD);
+    LinearAlgebra::TpetraWrappers::Vector<double> u;
+    u.reinit(complete_index_set(dim), MPI_COMM_WORLD);
+    f.trilinos_rcp()->putScalar(1.0);
+    A.compress(VectorOperation::insert);
+    f.compress(VectorOperation::insert);
+    u.compress(VectorOperation::insert);
+
+    GrowingVectorMemory<LinearAlgebra::TpetraWrappers::Vector<double>> mem;
+    SolverGMRES<LinearAlgebra::TpetraWrappers::Vector<double>> solver(control,
+                                                                      mem);
+    PreconditionIdentity                                       preconditioner;
+
+    deallog << "Solver type: " << typeid(solver).name() << std::endl;
+
+    check_solver_within_range(solver.solve(A, u, f, preconditioner),
+                              control.last_step(),
+                              74,
+                              76);
+  }
+}

tests/trilinos_tpetra/deal_solver_03.output

@@ -0,0 +1,4 @@
+
+DEAL::Size 32 Unknowns 961
+DEAL::Solver type: N6dealii11SolverGMRESINS_13LinearAlgebra14TpetraWrappers6VectorIdNS_11MemorySpace4HostEEEEE
+DEAL::Solver stopped within 74 - 76 iterations