Homec4science
Diffusion µSpectre (18_update_pfft) Browse

muSpectre/tests/header_test_tensor_algebra.cc
97c846c0c24518_update_pfft

Last Change
Raw File

header_test_tensor_algebra.cc
View Options

/**
* @file header_test_tensor_algebra.cc
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 05 Nov 2017
*
* @brief Tests for the tensor algebra functions
*
* Copyright © 2017 Till Junge
*
* µSpectre is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µSpectre is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include <iomanip>
#include <unsupported/Eigen/CXX11/Tensor>
#include "common/tensor_algebra.hh"
#include "tests.hh"
#include "tests/test_goodies.hh"
namespace muSpectre {
BOOST_AUTO_TEST_SUITE(tensor_algebra)
auto TerrNorm = [](auto && t){
return Eigen::Tensor<Real, 0>(t.abs().sum())();
};
/* ---------------------------------------------------------------------- */
BOOST_AUTO_TEST_CASE(tensor_outer_product_test) {
constexpr Dim_t dim{2};
Eigen::TensorFixedSize<Real, Eigen::Sizes<dim,dim>> A, B;
// use prime numbers so that every multiple is uniquely identifiable
A.setValues({{1, 2}, {3 ,7}});
B.setValues({{11, 13}, {17 ,19}});
Eigen::TensorFixedSize<Real, Eigen::Sizes<dim, dim, dim, dim>> Res1, Res2, Res3;
for (Dim_t i = 0; i < dim; ++i) {
for (Dim_t j = 0; j < dim; ++j) {
for (Dim_t k = 0; k < dim; ++k) {
for (Dim_t l = 0; l < dim; ++l) {
Res1(i, j, k, l) = A(i, j)*B(k, l);
Res2(i, j, k, l) = A(i, k)*B(j, l);
Res3(i, j, k, l) = A(i, l)*B(j, k);
}
}
}
}
Real error = TerrNorm(Res1 - Tensors::outer<dim>(A,B));
BOOST_CHECK_LT(error, tol);
error = TerrNorm(Res2 - Tensors::outer_under<dim>(A, B));
BOOST_CHECK_LT(error, tol);
error = TerrNorm(Res3 - Tensors::outer_over<dim>(A, B));
if (error > tol) {
std::cout << "reference:" << std::endl
<< Res3 << std::endl;
std::cout << "result:" << std::endl
<< Tensors::outer_over<dim>(A, B) << std::endl;
std::cout << "A:" << std::endl
<< A << std::endl;
std::cout << "B" << std::endl
<< B << std::endl;
decltype(Res3) tmp = Tensors::outer_over<dim>(A, B);
for (Dim_t i = 0; i < dim; ++i) {
for (Dim_t j = 0; j < dim; ++j) {
for (Dim_t k = 0; k < dim; ++k) {
for (Dim_t l = 0; l < dim; ++l) {
std::cout << "for (" << i << ", " << j << ", " << k << ", " << l
<< "), ref: " << std::setw(3)<< Res3(i,j,k,l)
<< ", res: " << std::setw(3)<< tmp(i,j,k,l)
<< std::endl;
}
}
}
}
}
BOOST_CHECK_LT(error, tol);
error = TerrNorm(Res3 - Tensors::outer<dim>(A, B));
BOOST_CHECK_GT(error, tol);
};
BOOST_FIXTURE_TEST_CASE_TEMPLATE(outer_products, Fix,
testGoodies::dimlist, Fix) {
constexpr auto dim{Fix::dim};
using T2 = Tensors::Tens2_t<dim>;
using M2 = Matrices::Tens2_t<dim>;
using Map2 = Eigen::Map<const M2>;
using T4 = Tensors::Tens4_t<dim>;
using M4 = Matrices::Tens4_t<dim>;
using Map4 = Eigen::Map<const M4>;
T2 A, B;
T4 RT;
A.setRandom();
B.setRandom();
Map2 Amap(A.data());
Map2 Bmap(B.data());
M2 C, D;
M4 RM;
C = Amap;
D = Bmap;
auto error =[](const T4& A, const M4& B) {return (B - Map4(A.data())).norm();};
// Check outer product
RT = Tensors::outer<dim>(A, B);
RM = Matrices::outer(C, D);
BOOST_CHECK_LT(error(RT, RM), tol);
// Check outer_under product
RT = Tensors::outer_under<dim>(A, B);
RM = Matrices::outer_under(C, D);
BOOST_CHECK_LT(error(RT, RM), tol);
// Check outer_over product
RT = Tensors::outer_over<dim>(A, B);
RM = Matrices::outer_over(C, D);
BOOST_CHECK_LT(error(RT, RM), tol);
}
BOOST_AUTO_TEST_CASE(tensor_multiplication) {
constexpr Dim_t dim{2};
using Strain_t = Eigen::TensorFixedSize<Real, Eigen::Sizes<dim, dim>>;
Strain_t A, B;
A.setValues({{1, 2}, {3 ,7}});
B.setValues({{11, 13}, {17 ,19}});
Strain_t FF1 = A*B; // element-wise multiplication
std::array<Eigen::IndexPair<int>, 1> prod_dims
{ Eigen::IndexPair<int>{1, 0}};
Strain_t FF2 = A.contract(B, prod_dims); // correct option 1
Strain_t FF3;
using Mat_t = Eigen::Map<Eigen::Matrix<Real, dim, dim>>;
// following only works for evaluated tensors (which already have data())
Mat_t(FF3.data()) = Mat_t(A.data()) * Mat_t(B.data());
Strain_t ref;
ref.setZero();
for (Dim_t i = 0; i < dim; ++i) {
for (Dim_t j = 0; j < dim; ++j) {
for (Dim_t a = 0; a < dim; ++a) {
ref(i,j) += A(i, a) * B(a, j);
}
}
}
using Strain_tw = Eigen::TensorFixedSize<Real, Eigen::Sizes<dim+1, dim+1>>;
Strain_tw C;
C.setConstant(100);
//static_assert(!std::is_convertible<Strain_t, Strain_tw>::value,
// "Tensors not size-protected");
if (std::is_convertible<Strain_t, Strain_tw>::value) {
//std::cout << "this is not good, should I abandon Tensors?";
}
// this test seems useless. I use to detect if Eigen changed the
// default tensor product
Real error = TerrNorm(FF1-ref);
if (error < tol) {
std::cout << "A =" << std::endl
<< A << std::endl;
std::cout << "B =" << std::endl
<< B << std::endl;
std::cout << "FF1 =" << std::endl
<< FF1 << std::endl;
std::cout << "ref =" << std::endl
<< ref << std::endl;
}
BOOST_CHECK_GT(error, tol);
error = TerrNorm(FF2-ref);
if (error > tol) {
std::cout << "A =" << std::endl
<< A << std::endl;
std::cout << "B =" << std::endl
<< B << std::endl;
std::cout << "FF2 =" << std::endl
<< FF2 << std::endl;
std::cout << "ref =" << std::endl
<< ref << std::endl;
}
BOOST_CHECK_LT(error, tol);
error = TerrNorm(FF3-ref);
if (error > tol) {
std::cout << "A =" << std::endl
<< A << std::endl;
std::cout << "B =" << std::endl
<< B << std::endl;
std::cout << "FF3 =" << std::endl
<< FF3 << std::endl;
std::cout << "ref =" << std::endl
<< ref << std::endl;
}
BOOST_CHECK_LT(error, tol);
}
BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tensmult, Fix,
testGoodies::dimlist, Fix) {
using namespace Matrices;
constexpr Dim_t dim{Fix::dim};
using T4 = Tens4_t<dim>;
using T2 = Tens2_t<dim>;
using V2 = Eigen::Matrix<Real, dim*dim, 1>;
T4 C;
C.setRandom();
T2 E;
E.setRandom();
Eigen::Map<const V2> Ev(E.data());
T2 R = tensmult(C, E);
auto error = (Eigen::Map<const V2>(R.data())-C*Ev).norm();
BOOST_CHECK_LT(error, tol);
}
BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tracer, Fix,
testGoodies::dimlist, Fix) {
using namespace Matrices;
constexpr Dim_t dim{Fix::dim};
using T2 = Tens2_t<dim>;
auto tracer = Itrac<dim>();
T2 F;
F.setRandom();
auto Ftrac = tensmult(tracer, F);
auto error = (Ftrac-F.trace()*F.Identity()).norm();
BOOST_CHECK_LT(error, tol);
}
BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_identity, Fix,
testGoodies::dimlist, Fix) {
using namespace Matrices;
constexpr Dim_t dim{Fix::dim};
using T2 = Tens2_t<dim>;
auto ident = Iiden<dim>();
T2 F;
F.setRandom();
auto Fiden = tensmult(ident, F);
auto error = (Fiden-F).norm();
BOOST_CHECK_LT(error, tol);
}
BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_transposer, Fix,
testGoodies::dimlist, Fix) {
using namespace Matrices;
constexpr Dim_t dim{Fix::dim};
using T2 = Tens2_t<dim>;
auto trnst = Itrns<dim>();
T2 F;
F.setRandom();
auto Ftrns = tensmult(trnst, F);
auto error = (Ftrns-F.transpose()).norm();
BOOST_CHECK_LT(error, tol);
}
BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_symmetriser, Fix,
testGoodies::dimlist, Fix) {
using namespace Matrices;
constexpr Dim_t dim{Fix::dim};
using T2 = Tens2_t<dim>;
auto symmt = Isymm<dim>();
T2 F;
F.setRandom();
auto Fsymm = tensmult(symmt, F);
auto error = (Fsymm-.5*(F+F.transpose())).norm();
BOOST_CHECK_LT(error, tol);
}
BOOST_AUTO_TEST_SUITE_END();
} // muSpectre
c4science · Help