numeric_vector_test.h

Go to the documentation of this file.

#ifndef NUMERIC_VECTOR_TEST_H
#define NUMERIC_VECTOR_TEST_H

// test includes
#include "test_comm.h"

// libMesh includes
#include <libmesh/parallel.h>
#include <libmesh/fuzzy_equals.h>
#include <libmesh/int_range.h>

#include "libmesh_cppunit.h"

#include <memory>

#define NUMERICVECTORTEST                       \
  CPPUNIT_TEST( testLocalize );                 \
  CPPUNIT_TEST( testLocalizeBase );             \
  CPPUNIT_TEST( testLocalizeIndices );          \
  CPPUNIT_TEST( testLocalizeIndicesBase );      \
  CPPUNIT_TEST( testLocalizeToOne );            \
  CPPUNIT_TEST( testLocalizeToOneBase );        \
  CPPUNIT_TEST( testNorms );                    \
  CPPUNIT_TEST( testNormsBase );                \
  CPPUNIT_TEST( testOperations );               \
  CPPUNIT_TEST( testOperationsBase );           \
  CPPUNIT_TEST( testWriteAndRead );             \
  CPPUNIT_TEST( testWriteAndReadBase );

// Tests to add manually, on subclasses supporting them
//  CPPUNIT_TEST( testSubvectors );
//  CPPUNIT_TEST( testSubvectorsBase );


template <class DerivedClass>
class NumericVectorTest : public CppUnit::TestCase {

protected:
  libMesh::Parallel::Communicator *my_comm = nullptr;

  std::string libmesh_suite_name;

  unsigned int block_size, local_size, global_size;

public:
  void setUp()
  {
    // By default we'll use the whole communicator in parallel;
    // Serial-only NumericVector subclasses will need to override
    // this and set something else first.
    if (!my_comm)
      my_comm = TestCommWorld;

    block_size  = 10;

    // a different size on each processor.
    local_size  = block_size +
      static_cast<unsigned int>(my_comm->rank());

    global_size = 0;
    for (libMesh::processor_id_type p=0; p<my_comm->size(); p++)
      global_size += (block_size + static_cast<unsigned int>(p));
  }

  void tearDown()
  {}

  template <class Base, class Derived>
  void Operations()
  {
    auto v_ptr = std::make_unique<Derived>(*my_comm, global_size, local_size);
    Base & v = *v_ptr;

    const libMesh::dof_id_type
      first = v.first_local_index(),
      last  = v.last_local_index();

    for (libMesh::dof_id_type n=first; n != last; n++)
      v.set (n, static_cast<libMesh::Number>(n+1));
    v.close();

    auto v_clone = v.clone();
    auto & vorig = *v_clone;

    CPPUNIT_ASSERT(relative_fuzzy_equals(v, vorig));

    v += v;

    CPPUNIT_ASSERT(!relative_fuzzy_equals(v, vorig));

    for (libMesh::dof_id_type n=first; n != last; n++)
      LIBMESH_ASSERT_NUMBERS_EQUAL
        (v(n), libMesh::Real(2*n+2),
         libMesh::TOLERANCE*libMesh::TOLERANCE);

    v *= v;

    for (libMesh::dof_id_type n=first; n != last; n++)
      LIBMESH_ASSERT_NUMBERS_EQUAL
        (v(n), libMesh::Real(4*n*n+8*n+4),
         libMesh::TOLERANCE*libMesh::TOLERANCE);

    v -= vorig;

    for (libMesh::dof_id_type n=first; n != last; n++)
      LIBMESH_ASSERT_NUMBERS_EQUAL
        (v(n), libMesh::Real(4*n*n+7*n+3),
         libMesh::TOLERANCE*libMesh::TOLERANCE);

    v /= vorig;

    for (libMesh::dof_id_type n=first; n != last; n++)
      LIBMESH_ASSERT_NUMBERS_EQUAL
        (v(n), libMesh::Real(4*n+3),
         libMesh::TOLERANCE*libMesh::TOLERANCE);

    v.add(1);
    for (libMesh::dof_id_type n=first; n != last; n++)
      LIBMESH_ASSERT_NUMBERS_EQUAL
        (v(n), libMesh::Real(4*n+4),
         libMesh::TOLERANCE*libMesh::TOLERANCE);

    v.add(2, vorig);
    for (libMesh::dof_id_type n=first; n != last; n++)
      LIBMESH_ASSERT_NUMBERS_EQUAL
        (v(n), libMesh::Real(6*n+6),
         libMesh::TOLERANCE*libMesh::TOLERANCE);

    v.scale(1/libMesh::Real(3));
    for (libMesh::dof_id_type n=first; n != last; n++)
      LIBMESH_ASSERT_NUMBERS_EQUAL
        (v(n), libMesh::Real(2*n+2),
         libMesh::TOLERANCE*libMesh::TOLERANCE);

    v = 4;
    for (libMesh::dof_id_type n=first; n != last; n++)
      LIBMESH_ASSERT_NUMBERS_EQUAL
        (v(n), libMesh::Real(4),
         libMesh::TOLERANCE*libMesh::TOLERANCE);

    v = vorig;
    for (libMesh::dof_id_type n=first; n != last; n++)
      LIBMESH_ASSERT_NUMBERS_EQUAL
        (v(n), libMesh::Real(n+1),
         libMesh::TOLERANCE*libMesh::TOLERANCE);

    v.reciprocal();
    for (libMesh::dof_id_type n=first; n != last; n++)
      LIBMESH_ASSERT_NUMBERS_EQUAL
        (v(n), 1/libMesh::Real(n+1),
         libMesh::TOLERANCE*libMesh::TOLERANCE);

    LIBMESH_ASSERT_NUMBERS_EQUAL
      (v.dot(vorig), libMesh::Real(global_size),
       libMesh::TOLERANCE*libMesh::TOLERANCE);
  }


  template <class Base, class Derived>
  void Subvectors()
  {
    auto v_ptr = std::make_unique<Derived>(*my_comm, global_size, local_size);
    Base & v = *v_ptr;

    const libMesh::dof_id_type
      first = v.first_local_index(),
      last  = v.last_local_index();

    const unsigned int sub_local_size = block_size/2;
    const unsigned int sub_global_size = sub_local_size * my_comm->size();

    for (auto supplying_global_rows : {false, true})
      for (auto get_not_create : {false, true})
        {
          // We don't have such an option?
          if (get_not_create && supplying_global_rows)
            continue;

          for (libMesh::dof_id_type n=first; n != last; n++)
            v.set (n, static_cast<libMesh::Number>(n+1));
          v.close();

          std::unique_ptr<Base> s_ptr =
            std::make_unique<Derived>(*my_comm, sub_global_size,
                                      sub_local_size);

          const libMesh::dof_id_type
            sub_first = s_ptr->first_local_index(),
            sub_last  = s_ptr->last_local_index();

          CPPUNIT_ASSERT_EQUAL
            (sub_first, libMesh::dof_id_type(sub_local_size *
                                             my_comm->rank()));

          CPPUNIT_ASSERT_EQUAL
            (sub_last, libMesh::dof_id_type(sub_local_size *
                                            (my_comm->rank()+1)));

          std::vector<libMesh::dof_id_type> rows;
          for (auto i : libMesh::make_range(sub_local_size))
            rows.push_back(first+2*i);

          if (supplying_global_rows)
            my_comm->allgather(rows);

          if (get_not_create)
            {
              libMesh::NumericVector<libMesh::Number> * subvec_ptr =
                v.get_subvector(rows).release();
              s_ptr.reset(dynamic_cast<Base *>(subvec_ptr));
            }
          else
            v.create_subvector(*s_ptr, rows, supplying_global_rows);

          Base & s = *s_ptr;

          for (auto i : libMesh::make_range(sub_local_size))
            LIBMESH_ASSERT_NUMBERS_EQUAL
              (s(sub_first+i), libMesh::Real(first+2*i+1),
               libMesh::TOLERANCE*libMesh::TOLERANCE);

          for (libMesh::dof_id_type n=sub_first; n != sub_last; n++)
            s.set(n, (my_comm->rank()+1) * n);
          s.close();

          // If we "get" a subvector then we're set up to "restore" it
          // afterwards; if not then we're done.
          if (!get_not_create)
            continue;

          v.restore_subvector(std::move(s_ptr), rows);

          for (libMesh::dof_id_type n=first; n != last; n++)
            {
              const libMesh::dof_id_type offset = n - first;
              if (offset < 2*sub_local_size && !(offset%2))
                LIBMESH_ASSERT_NUMBERS_EQUAL
                  (v(n),
                   libMesh::Real((my_comm->rank()+1) *
                                 ((offset/2) + sub_first)),
                   libMesh::TOLERANCE*libMesh::TOLERANCE);
              else
                LIBMESH_ASSERT_NUMBERS_EQUAL
                  (v(n), libMesh::Real(n+1),
                   libMesh::TOLERANCE*libMesh::TOLERANCE);
            }
        }
  }

  template <class Base, class Derived>
  void WriteAndRead()
  {
    auto v_ptr = std::make_unique<Derived>(*my_comm, global_size, local_size);
    Base & v = *v_ptr;

    const libMesh::dof_id_type
      first = v.first_local_index(),
      last  = v.last_local_index();

    for (libMesh::dof_id_type n=first; n != last; n++)
      v.set (n, static_cast<libMesh::Number>(n+1));
    v.close();

    v.print_matlab(libmesh_suite_name+"_vector.m");

    // Let's make sure we don't have any race conditions; we have
    // multiple SparseMatrix subclasses that might be trying to read
    // and write the same file.

    TestCommWorld->barrier();

#ifdef LIBMESH_USE_COMPLEX_NUMBERS
    // We're not supporting complex reads quite yet
    return;
#endif

    auto v_new_ptr = std::make_unique<Derived>(*my_comm, global_size, local_size);
    Base & v_new = *v_new_ptr;

    v_new.read_matlab(libmesh_suite_name+"_vector.m");

    TestCommWorld->barrier();

    CPPUNIT_ASSERT_EQUAL(v.l1_norm_diff(v_new), libMesh::Real(0));
  }

  template <class Base, class Derived>
  void Norms()
  {
    auto v_ptr = std::make_unique<Derived>(*my_comm, global_size, local_size);
    Base & v = *v_ptr;

    const libMesh::dof_id_type
      first = v.first_local_index(),
      last  = v.last_local_index();

    for (libMesh::dof_id_type n=first; n != last; n++)
      v.set (n, -static_cast<libMesh::Number>(n));
    v.close();
    for (libMesh::dof_id_type n=first; n != last; n++)
      v.add (n, -static_cast<libMesh::Number>(n));
    v.close();

    const libMesh::Real exact_l1 =
      global_size * (global_size-libMesh::Real(1));
    const libMesh::Real exact_l2 =
      std::sqrt(libMesh::Real(global_size-1) *
                (2*global_size) *
                (2*global_size-1) / 3);
    LIBMESH_ASSERT_NUMBERS_EQUAL(v.sum(), -exact_l1,
                            libMesh::TOLERANCE*libMesh::TOLERANCE);
    LIBMESH_ASSERT_FP_EQUAL(v.l1_norm(), exact_l1,
                            libMesh::TOLERANCE*libMesh::TOLERANCE);
    LIBMESH_ASSERT_FP_EQUAL(v.l2_norm(), exact_l2,
                            libMesh::TOLERANCE*libMesh::TOLERANCE);
    LIBMESH_ASSERT_FP_EQUAL(v.linfty_norm(), 2*(global_size-libMesh::Real(1)),
                            libMesh::TOLERANCE*libMesh::TOLERANCE);

    auto u_ptr = std::make_unique<Derived>(*my_comm, global_size, local_size);
    Base & u = *u_ptr;
    for (libMesh::dof_id_type n=first; n != last; n++)
      u.set (n, -static_cast<libMesh::Number>(n * n));
    u.close();

    Derived diff_derived(*my_comm, global_size, local_size);
    Base & diff = diff_derived;
    diff = u;
    diff -= v;
    diff.close();

    // Use a relative tolerance here; the norms are O(1e5) when our
    // processor count gets big enough, at which point O(1e-12) is
    // testing for exact equality...
    const auto diff_norm = diff.l2_norm();
    const auto norm_diff = u.l2_norm_diff(v);
    LIBMESH_ASSERT_FP_EQUAL(diff_norm, norm_diff,
                            (std::abs(diff_norm)+std::abs(norm_diff)) *
                            libMesh::TOLERANCE*libMesh::TOLERANCE);
  }

  template <class Base, class Derived>
  void Localize(bool to_one=false)
  {
    const libMesh::processor_id_type root_pid = 0;

    {
      auto v_ptr = std::make_unique<Derived>(*my_comm, global_size, local_size);
      Base & v = *v_ptr;
      std::vector<libMesh::Number> l(global_size);

      const libMesh::dof_id_type
        first = v.first_local_index(),
        last  = v.last_local_index();

      for (libMesh::dof_id_type n=first; n != last; n++)
        v.set (n, static_cast<libMesh::Number>(n));
      v.close();
      for (libMesh::dof_id_type n=first; n != last; n++)
        v.add (n, static_cast<libMesh::Number>(n));
      v.close();

      if (!to_one)
        v.localize(l);
      else
        v.localize_to_one(l,root_pid);

      if (!to_one || my_comm->rank() == root_pid)
        // Yes I really mean v.size()
        for (libMesh::dof_id_type i=0; i<v.size(); i++)
          LIBMESH_ASSERT_NUMBERS_EQUAL
            (2.*i, l[i], libMesh::TOLERANCE*libMesh::TOLERANCE);

      for (libMesh::dof_id_type n=first; n != last; n++)
      {
        const auto value = static_cast<libMesh::Number>(n);
        v.insert (&value, std::vector<libMesh::numeric_index_type>({n}));
      }
      v.close();

      if (!to_one)
        v.localize(l);
      else
        v.localize_to_one(l,root_pid);

      if (!to_one || my_comm->rank() == root_pid)
        // Yes I really mean v.size()
        for (libMesh::dof_id_type i=0; i<v.size(); i++)
          LIBMESH_ASSERT_NUMBERS_EQUAL
            (libMesh::Real(i), l[i], libMesh::TOLERANCE*libMesh::TOLERANCE);
    }
  }


  template <class Base, class Derived>
  void LocalizeIndices()
  {
    {
      auto v_ptr = std::make_unique<Derived>(*my_comm, global_size, local_size);
      Base & v = *v_ptr;

      // Let's try pulling the same number of entries from each processor
      std::vector<libMesh::Number> values(block_size * my_comm->size());
      std::vector<libMesh::dof_id_type> indices;
      indices.reserve(block_size * my_comm->size());

      const libMesh::dof_id_type
        first = v.first_local_index(),
        last  = v.last_local_index();

      for (libMesh::dof_id_type n=first; n != last; n++)
        v.set (n, static_cast<libMesh::Number>(n));
      v.close();

      libMesh::dof_id_type end_index = 0;
      for (libMesh::processor_id_type p=0; p<my_comm->size(); p++)
        {
          end_index += block_size + p;
          for (unsigned int j = 0; j != block_size; ++j)
            indices.push_back(end_index-j-1);
        }

      v.localize(values, indices);

      end_index = 0;
      for (libMesh::processor_id_type p=0; p<my_comm->size(); p++)
        {
          end_index += block_size + p;
          for (unsigned int j = 0; j != block_size; ++j)
            LIBMESH_ASSERT_NUMBERS_EQUAL
              (values[p*block_size+j], libMesh::Real(end_index-j-1),
               libMesh::TOLERANCE*libMesh::TOLERANCE);
        }
    }
  }

  void testLocalize()
  {
    LOG_UNIT_TEST;

    Localize<DerivedClass,DerivedClass>();
  }

  void testLocalizeBase()
  {
    LOG_UNIT_TEST;

    Localize<libMesh::NumericVector<libMesh::Number>,DerivedClass>();
  }

  void testLocalizeToOne()
  {
    LOG_UNIT_TEST;

    Localize<DerivedClass,DerivedClass >(true);
  }

  void testLocalizeToOneBase()
  {
    LOG_UNIT_TEST;

    Localize<libMesh::NumericVector<libMesh::Number>,DerivedClass>(true);
  }

  void testLocalizeIndices()
  {
    LOG_UNIT_TEST;

    LocalizeIndices<DerivedClass,DerivedClass >();
  }

  void testLocalizeIndicesBase()
  {
    LOG_UNIT_TEST;

    LocalizeIndices<libMesh::NumericVector<libMesh::Number>,DerivedClass>();
  }

  void testNorms()
  {
    LOG_UNIT_TEST;

    Norms<DerivedClass,DerivedClass >();
  }

  void testNormsBase()
  {
    LOG_UNIT_TEST;

    Norms<libMesh::NumericVector<libMesh::Number>,DerivedClass>();
  }

  void testOperations()
  {
    LOG_UNIT_TEST;

    Operations<DerivedClass,DerivedClass >();
  }

  void testOperationsBase()
  {
    LOG_UNIT_TEST;

    Operations<libMesh::NumericVector<libMesh::Number>,DerivedClass>();
  }

  void testSubvectors()
  {
    LOG_UNIT_TEST;

    Subvectors<DerivedClass,DerivedClass >();
  }

  void testSubvectorsBase()
  {
    LOG_UNIT_TEST;

    Subvectors<libMesh::NumericVector<libMesh::Number>,DerivedClass>();
  }

  void testWriteAndRead()
  {
    LOG_UNIT_TEST;

    WriteAndRead<DerivedClass,DerivedClass >();
  }

  void testWriteAndReadBase()
  {
    LOG_UNIT_TEST;

    WriteAndRead<libMesh::NumericVector<libMesh::Number>,DerivedClass>();
  }
};

#endif