dbeurle / neon

A finite element code
Other
12 stars 8 forks source link

Voigt vs Mandel (or Kelvin) notation #142

Open dbeurle opened 5 years ago

dbeurle commented 5 years ago

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));
                           });
  }
shadisharba commented 5 years ago

I ran the same test and got some interesting results. The voigt_contraction requires only two thirds the time to run tensor_contraction.

2019-01-20 17:19:02
Running ./a.out
Run on (8 X 4000 MHz CPU s)
CPU Caches:
  L1 Data 32K (x4)
  L1 Instruction 32K (x4)
  L2 Unified 256K (x4)
  L3 Unified 8192K (x1)
Load Average: 0.85, 1.13, 1.08
------------------------------------------------------------------------------
Benchmark                   Time             CPU   Iterations UserCounters...
------------------------------------------------------------------------------
tensor_contraction      29366 ns        29314 ns        21955 norms=21.955k
voigt_contraction       20264 ns        20232 ns        34895 norms=34.895k
dbeurle commented 5 years ago

Okay your results make more sense. I'll test on my workstation as well when I get home.

augustobopsinborges commented 3 years ago

Are you considering that Mandel notation consists in representing 2nd order tensors as 3×3 matrices? As far as I know, Mandel notation expresses also 2nd order tensors in vector form, but in the following form: T{ij} = T{ji} → [T{11} T{22} T{33} \sqrt{2}T{23} \sqrt{2}T{13} \sqrt{2}T{12}]^t

You seem not clear what you are considering Mandel notation. Can you clarify what are you considering Mandel notation?