Skip to content

Voigt vs Mandel (or Kelvin) notation #142

@dbeurle

Description

@dbeurle

We current store second order tensors as 3x3 matrices, regardless of symmetric considerations. Some memory (or theoretically computation) reduction can be achieved from using Voigt or Mandel notation.

Some data on a conversion:

Benchmark Number of norm operations
Voigt (no 2 correction) 114787ns
Tensor 103020ns
Voigt (with 2 correction) 201364ns

Benchmark code

Running ./run_benchmark
Run on (4 X 3500 MHz CPU s)
CPU Caches:
  L1 Data 32K (x2)
  L1 Instruction 32K (x2)
  L2 Unified 256K (x2)
  L3 Unified 3072K (x1)
***WARNING*** CPU scaling is enabled, the benchmark real time measurements may be noisy and will incur extra overhead.
--------------------------------------------------------------------------
Benchmark                   Time           CPU Iterations UserCounters...
--------------------------------------------------------------------------
tensor_contraction     115418 ns     114787 ns       5990 norms=5.99k
voigt_contraction      103547 ns     103020 ns       6670 norms=6.67k

Provided by the code:

#include <benchmark/benchmark.h>

#include <cmath>
#include <numeric>

#include <eigen3/Eigen/Dense>

using vector6 = Eigen::Matrix<double, 6, 1>;
using matrix3 = Eigen::Matrix<double, 3, 3>;

static void tensor_contraction(benchmark::State &state) {

  std::vector<matrix3> as(10'000), bs(10'000);

  std::generate(begin(as), end(as), []() { return matrix3::Random(3, 3); });
  std::generate(begin(bs), end(bs), []() { return matrix3::Random(3, 3); });

  for (auto _ : state) {
    [[maybe_unused]] volatile double result =
        std::inner_product(begin(as), end(as), begin(bs), 0.0, std::plus<>(),
                           [](auto const &a, auto const &b) {
                             return std::sqrt((a.array() * b.array()).sum());
                           });
  }
  state.counters["norms"] = state.iterations();
}

static void voigt_contraction(benchmark::State &state) {

  std::vector<vector6> as(10'000), bs(10'000);

  std::generate(begin(as), end(as), []() { return vector6::Random(6, 1); });
  std::generate(begin(bs), end(bs), []() { return vector6::Random(6, 1); });

  for (auto _ : state) {
    [[maybe_unused]] volatile double result = std::inner_product(
        begin(as), end(as), begin(bs), 0.0, std::plus<>(),
        [](auto const &a, auto const &b) { return std::sqrt(a.dot(b)); });
  }

  state.counters["norms"] = state.iterations();
}

BENCHMARK(tensor_contraction);
BENCHMARK(voigt_contraction);

BENCHMARK_MAIN();

and

  for (auto _ : state) {
    [[maybe_unused]] volatile double result =
        std::inner_product(begin(as), end(as), begin(bs), 0.0, std::plus<>(),
                           [](vector6 &a, vector6 &b) {
                             a.tail<3>() *= 2.0;
                             b.tail<3>() *= 2.0;

                             return std::sqrt(a.dot(b));
                           });
  }

Metadata

Metadata

Assignees

No one assigned

    Labels

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions