Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/tests/test_integration.cpp
acb92043a41dfeatures/netcdf-mpi

Last Change
Raw File

test_integration.cpp
View Options

/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "fftransform.hh"
#include "integration/accumulator.hh"
#include "test.hh"
#include <cmath>
#include <expolit/expolit>
/* -------------------------------------------------------------------------- */
using namespace tamaas;
namespace ex = expolit;
static Real tol = 1e-13;
class AccumulatorTest : public ::testing::Test {
protected:
void SetUp() override {
nodal_values.resize(n);
for (auto&& grid : nodal_values) {
grid.setNbComponents(6);
grid.resize({nh, nh});
grid = 1;
}
accumulator.makeUniformMesh(nodal_values.size(), size);
wavevectors = FFTransform<Real, 2>::computeFrequencies<true>({nh, nh});
}
UInt n = 10, nh = 10;
Real size = 3;
std::vector<GridHermitian<Real, 2>> nodal_values;
Grid<Real, 2> wavevectors;
Accumulator<model_type::volume_2d, SymMatrixProxy<Complex, 3>> accumulator;
};
TEST_F(AccumulatorTest, uniformMesh) {
auto& pos = this->accumulator.nodePositions();
std::vector<Real> sol(pos.size());
std::iota(sol.begin(), sol.end(), 0);
std::transform(sol.begin(), sol.end(), sol.begin(),
[&](auto& x) { return size * x / (sol.size() - 1); });
ASSERT_TRUE(compare(sol, pos, AreFloatEqual())) << "Uniform mesh fail";
}
// Litteral type
using Q = ex::Litteral<struct Q_>;
TEST_F(AccumulatorTest, forward) {
// Setup for layer checking
std::vector<int> layers_seen, layers_to_see(n);
std::iota(layers_to_see.begin(), layers_to_see.end(), 0);
// Setup for integral checking
constexpr Q q;
constexpr auto z = ex::Polynomial<Real, 1>({0, 1});
constexpr auto g0 = ex::exp(q * z);
constexpr auto g1 = q * z * ex::exp(q * z);
for (auto&& tuple :
accumulator.forward(this->nodal_values, this->wavevectors)) {
auto&& l = std::get<0>(tuple);
auto&& xl_ = std::get<1>(tuple);
auto&& acc_g0 = std::get<2>(tuple);
auto&& acc_g1 = std::get<3>(tuple);
layers_seen.push_back(l);
// fundamental mode
auto testing = acc_g0(0, 0, 0);
ASSERT_NEAR(testing.real(), xl_, tol) << "Fundamental mode G0 fail";
testing = acc_g1(0, 0, 0);
ASSERT_NEAR(testing.real(), 0, tol) << "Fundamental mode G1 fail";
// other modes
testing = acc_g0(1, 0, 0);
ASSERT_NEAR(testing.real(), acc_g0(0, 1, 0).real(), tol) << "Symmetry fail";
ASSERT_NEAR(testing.imag(), acc_g0(0, 1, 0).imag(), tol) << "Symmetry fail";
testing = acc_g0(1, 1, 0);
Real sqrt_two = std::sqrt(Real(2));
auto ref_value = ex::definite_integral(
std::pair<Real, Real>(0, xl_),
ex::substitute<Q::tag>(ex::Constant<Real>({sqrt_two}), g0));
ASSERT_NEAR(testing.real(), ref_value, tol)
<< "Integration of exp(qy)*ϕ(y) fail";
testing = acc_g1(1, 1, 0);
ref_value = ex::definite_integral(
std::pair<Real, Real>(0, xl_),
ex::substitute<Q::tag>(ex::Constant<Real>({sqrt_two}), g1));
ASSERT_NEAR(testing.real(), ref_value, tol)
<< "Integration of qy*exp(qy)*ϕ(y) fail";
}
ASSERT_TRUE(compare(layers_seen, layers_to_see))
<< "Did not see correct layers";
}
TEST_F(AccumulatorTest, backward) {
// Setup for layer checking
std::vector<int> layers_seen, layers_to_see(n);
std::iota(layers_to_see.rbegin(), layers_to_see.rend(), 0);
// Setup for integral checking
constexpr Q q;
constexpr auto z = ex::Polynomial<Real, 1>({0, 1});
constexpr auto g0 = ex::exp(q * (z * (-1)));
constexpr auto g1 = q * z * ex::exp(q * ((-1) * z));
for (auto&& tuple : accumulator.backward(nodal_values, wavevectors)) {
auto&& l = std::get<0>(tuple);
auto&& xl_ = std::get<1>(tuple);
auto&& acc_g0 = std::get<2>(tuple);
auto&& acc_g1 = std::get<3>(tuple);
layers_seen.push_back(l);
// fundamental mode
auto testing = acc_g0(0, 0, 0);
ASSERT_NEAR(testing.real(), size - xl_, tol) << "Fundamental mode G0 fail";
testing = acc_g1(0, 0, 0);
ASSERT_NEAR(testing.real(), 0, tol) << "Fundamental mode G1 fail";
// other modes
testing = acc_g0(1, 0, 0);
ASSERT_NEAR(testing.real(), acc_g0(0, 1, 0).real(), tol) << "Symmetry fail";
ASSERT_NEAR(testing.imag(), acc_g0(0, 1, 0).imag(), tol) << "Symmetry fail";
testing = acc_g0(1, 1, 0);
auto ref_value = ex::definite_integral(
std::make_pair(xl_, size),
ex::substitute<Q::tag>(ex::Constant<Real>({std::sqrt(2)}), g0));
ASSERT_NEAR(testing.real(), ref_value, tol)
<< "Integration of exp(-qy)*ϕ(y) fail";
testing = acc_g1(1, 1, 0);
ref_value = ex::definite_integral(
std::make_pair(xl_, size),
ex::substitute<Q::tag>(ex::Constant<Real>({std::sqrt(2)}), g1));
ASSERT_NEAR(testing.real(), ref_value, tol)
<< "Integration of qy*exp(-qy)*ϕ(y) fail";
}
ASSERT_TRUE(compare(layers_seen, layers_to_see))
<< "Did not see correct layers";
}
c4science · Help