Select Git revision
FPEManager.cpp
PolynomialReconstruction.cpp 20.76 KiB
#include <scheme/PolynomialReconstruction.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
#include <scheme/DiscreteFunctionVariant.hpp>
class MutableDiscreteFunctionDPkVariant
{
public:
using Variant = std::variant<DiscreteFunctionDPk<1, double>,
DiscreteFunctionDPk<1, TinyVector<1>>,
DiscreteFunctionDPk<1, TinyVector<2>>,
DiscreteFunctionDPk<1, TinyVector<3>>,
DiscreteFunctionDPk<1, TinyMatrix<1>>,
DiscreteFunctionDPk<1, TinyMatrix<2>>,
DiscreteFunctionDPk<1, TinyMatrix<3>>,
DiscreteFunctionDPk<2, double>,
DiscreteFunctionDPk<2, TinyVector<1>>,
DiscreteFunctionDPk<2, TinyVector<2>>,
DiscreteFunctionDPk<2, TinyVector<3>>,
DiscreteFunctionDPk<2, TinyMatrix<1>>,
DiscreteFunctionDPk<2, TinyMatrix<2>>,
DiscreteFunctionDPk<2, TinyMatrix<3>>,
DiscreteFunctionDPk<3, double>,
DiscreteFunctionDPk<3, TinyVector<1>>,
DiscreteFunctionDPk<3, TinyVector<2>>,
DiscreteFunctionDPk<3, TinyVector<3>>,
DiscreteFunctionDPk<3, TinyMatrix<1>>,
DiscreteFunctionDPk<3, TinyMatrix<2>>,
DiscreteFunctionDPk<3, TinyMatrix<3>>>;
#warning add DiscreteFunctionDPkVector to the variant
private:
Variant m_mutable_discrete_function_dpk;
public:
PUGS_INLINE
const Variant&
mutableDiscreteFunctionDPk() const
{
return m_mutable_discrete_function_dpk;
}
template <typename DiscreteFunctionDPkT>
PUGS_INLINE auto&&
get() const
{
static_assert(is_discrete_function_dpk_v<DiscreteFunctionDPkT>, "invalid template argument");
#ifndef NDEBUG
if (not std::holds_alternative<DiscreteFunctionDPkT>(this->m_mutable_discrete_function_dpk)) {
std::ostringstream error_msg;
error_msg << "invalid discrete function type\n";
error_msg << "- required " << rang::fgB::red << demangle<DiscreteFunctionDPkT>() << rang::fg::reset << '\n';
error_msg << "- contains " << rang::fgB::yellow
<< std::visit([](auto&& f) -> std::string { return demangle<decltype(f)>(); },
this->m_mutable_discrete_function_dpk)
<< rang::fg::reset;
throw NormalError(error_msg.str());
}
#endif // NDEBUG
return std::get<DiscreteFunctionDPkT>(this->mutableDiscreteFunctionDPk());
}
template <size_t Dimension, typename DataType>
MutableDiscreteFunctionDPkVariant(const DiscreteFunctionDPk<Dimension, DataType>& discrete_function_dpk)
: m_mutable_discrete_function_dpk{DiscreteFunctionDPk<Dimension, DataType>{discrete_function_dpk}} {
static_assert(std::is_same_v<std::remove_const_t<DataType>, double> or //
std::is_same_v<std::remove_const_t<DataType>, TinyVector<1, double>> or //
std::is_same_v<std::remove_const_t<DataType>, TinyVector<2, double>> or //
std::is_same_v<std::remove_const_t<DataType>, TinyVector<3, double>> or //
std::is_same_v<std::remove_const_t<DataType>, TinyMatrix<1, 1, double>> or //
std::is_same_v<std::remove_const_t<DataType>, TinyMatrix<2, 2, double>> or //
std::is_same_v<std::remove_const_t<DataType>, TinyMatrix<3, 3, double>>,
"DiscreteFunctionDPk with this DataType is not allowed in variant");
}
MutableDiscreteFunctionDPkVariant& operator=(MutableDiscreteFunctionDPkVariant&&) = default;
MutableDiscreteFunctionDPkVariant& operator=(const MutableDiscreteFunctionDPkVariant&) = default;
MutableDiscreteFunctionDPkVariant(const MutableDiscreteFunctionDPkVariant&) = default;
MutableDiscreteFunctionDPkVariant(MutableDiscreteFunctionDPkVariant&&) = default;
MutableDiscreteFunctionDPkVariant() = delete;
~MutableDiscreteFunctionDPkVariant() = default;
};
template <MeshConcept MeshType>
[[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
PolynomialReconstruction::_build(
const std::shared_ptr<const MeshType>& p_mesh,
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
{
const MeshType& mesh = *p_mesh;
using Rd = TinyVector<MeshType::Dimension>;
const size_t number_of_columns = [&]() {
size_t n = 0;
for (auto discrete_function_variant_p : discrete_function_variant_list) {
n += std::visit(
[](auto&& discrete_function) -> size_t {
using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
using data_type = std::decay_t<typename DiscreteFunctionT::data_type>;
if constexpr (std::is_same_v<data_type, double>) {
return 1;
} else if constexpr (is_tiny_vector_v<data_type> or is_tiny_matrix_v<data_type>) {
return data_type::Dimension;
} else {
throw UnexpectedError("unexpected data type " + demangle<data_type>());
}
} else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
return discrete_function.size();
} else {
throw UnexpectedError("unexpected discrete function type");
}
},
discrete_function_variant_p->discreteFunction());
}
return n;
}();
const size_t degree = 1;
const Stencil& stencil = StencilManager::instance().getStencil(mesh.connectivity());
auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
auto cell_is_owned = mesh.connectivity().cellIsOwned();
const size_t max_stencil_size = [&]() {
size_t max_size = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
auto stencil_cell_list = stencil[cell_id];
if (cell_is_owned[cell_id] and stencil_cell_list.size() > max_size) {
max_size = stencil_cell_list.size();
} }
return max_size;
}();
Kokkos::Experimental::UniqueToken<Kokkos::DefaultExecutionSpace::execution_space,
Kokkos::Experimental::UniqueTokenScope::Global>
tokens;
std::vector<MutableDiscreteFunctionDPkVariant> mutable_discrete_function_dpk_variant_list;
for (size_t i_discrete_function_variant = 0; i_discrete_function_variant < discrete_function_variant_list.size();
++i_discrete_function_variant) {
auto discrete_function_variant = discrete_function_variant_list[i_discrete_function_variant];
std::visit(
[&](auto&& discrete_function) {
using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
using DataType = std::remove_const_t<std::decay_t<typename DiscreteFunctionT::data_type>>;
mutable_discrete_function_dpk_variant_list.push_back(
DiscreteFunctionDPk<MeshType::Dimension, DataType>(p_mesh, degree));
} else {
throw NotImplementedError("discrete function type");
}
},
discrete_function_variant->discreteFunction());
}
SmallArray<SmallMatrix<double>> A_pool(Kokkos::DefaultExecutionSpace::concurrency());
SmallArray<SmallMatrix<double>> B_pool(Kokkos::DefaultExecutionSpace::concurrency());
SmallArray<SmallMatrix<double>> X_pool(Kokkos::DefaultExecutionSpace::concurrency());
for (size_t i = 0; i < A_pool.size(); ++i) {
A_pool[i] = SmallMatrix<double>(max_stencil_size, MeshType::Dimension);
B_pool[i] = SmallMatrix<double>(max_stencil_size, number_of_columns);
X_pool[i] = SmallMatrix<double>(MeshType::Dimension, number_of_columns);
}
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_j_id) {
if (cell_is_owned[cell_j_id]) {
const int32_t t = tokens.acquire();
auto stencil_cell_list = stencil[cell_j_id];
ShrinkMatrixView B(B_pool[t], stencil_cell_list.size());
for (size_t i_discrete_function_variant = 0;
i_discrete_function_variant < discrete_function_variant_list.size(); ++i_discrete_function_variant) {
const auto& discrete_function_variant = discrete_function_variant_list[i_discrete_function_variant];
std::visit(
[&](auto&& discrete_function) {
using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
size_t column_begin = 0;
const DataType& pj = discrete_function[cell_j_id];
for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
const CellId cell_i_id = stencil_cell_list[i];
const DataType& pi_pj = discrete_function[cell_i_id] - pj;
if constexpr (std::is_arithmetic_v<DataType>) {
B(i, column_begin) = pi_pj;
} else if constexpr (is_tiny_vector_v<DataType>) {
for (size_t kB = column_begin, k = 0; k < DataType::Dimension; ++k, ++kB) {
B(i, kB) = pi_pj[k];
}
} else if constexpr (is_tiny_matrix_v<DataType>) {
for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
B(i, column_begin + k * DataType::NumberOfColumns + l) = pi_pj(k, l); }
}
}
}
if constexpr (std::is_arithmetic_v<DataType>) {
++column_begin;
} else if constexpr (is_tiny_vector_v<DataType> or is_tiny_matrix_v<DataType>) {
column_begin += DataType::Dimension;
}
}
},
discrete_function_variant->discreteFunction());
}
// for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
// const CellId cell_i_id = stencil_cell_list[i];
// for (size_t j = 0; j < number_of_columns; ++j) {
// B(i, j) = values[cell_i_id][j] - values[cell_j_id][j];
// }
// }
ShrinkMatrixView A(A_pool[t], stencil_cell_list.size());
const Rd& Xj = xj[cell_j_id];
for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
const CellId cell_i_id = stencil_cell_list[i];
const Rd Xi_Xj = xj[cell_i_id] - Xj;
for (size_t l = 0; l < MeshType::Dimension; ++l) {
A(i, l) = Xi_Xj[l];
}
}
const SmallMatrix<double>& X = X_pool[t];
Givens::solveCollectionInPlace(A, X, B);
size_t column_begin = 0;
for (size_t i_dpk_variant = 0; i_dpk_variant < mutable_discrete_function_dpk_variant_list.size();
++i_dpk_variant) {
const auto& dpk_variant = mutable_discrete_function_dpk_variant_list[i_dpk_variant];
const auto& discrete_function_variant = discrete_function_variant_list[i_dpk_variant];
std::visit(
[&](auto&& dpk_function, auto&& p0_function) {
using DPkFunctionT = std::decay_t<decltype(dpk_function)>;
using P0FunctionT = std::decay_t<decltype(p0_function)>;
using DataType = std::remove_const_t<std::decay_t<typename DPkFunctionT::data_type>>;
using P0DataType = std::remove_const_t<std::decay_t<typename P0FunctionT::data_type>>;
if constexpr (std::is_same_v<DataType, P0DataType>) {
if constexpr (is_discrete_function_P0_v<P0FunctionT>) {
auto dpk_j = dpk_function.coefficients(cell_j_id);
dpk_j[0] = p0_function[cell_j_id];
if constexpr (std::is_arithmetic_v<DataType>) {
for (size_t i = 0; i < MeshType::Dimension; ++i) {
auto& dpk_j_ip1 = dpk_j[i + 1];
dpk_j_ip1 = X(i, column_begin);
}
++column_begin;
} else if constexpr (is_tiny_vector_v<DataType>) {
for (size_t i = 0; i < MeshType::Dimension; ++i) {
auto& dpk_j_ip1 = dpk_j[i + 1];
for (size_t k = 0; k < DataType::Dimension; ++k) {
dpk_j_ip1[k] = X(i, column_begin + k);
}
}
column_begin += DataType::Dimension; } else if constexpr (is_tiny_matrix_v<DataType>) {
for (size_t i = 0; i < MeshType::Dimension; ++i) {
auto& dpk_j_ip1 = dpk_j[i + 1];
for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
dpk_j_ip1(k, l) = X(i, column_begin + k * DataType::NumberOfColumns + l);
}
}
}
column_begin += DataType::Dimension;
} else {
throw NotImplementedError("unexpected data type");
}
} else {
throw NotImplementedError("non scalar P0 functions");
}
} else {
throw UnexpectedError("incompatible data types");
}
},
dpk_variant.mutableDiscreteFunctionDPk(), discrete_function_variant->discreteFunction());
}
tokens.release(t);
}
});
return {};
// parallel_for(
// mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_j_id) {
// if (cell_is_owned[cell_j_id]) {
// const int32_t t = tokens.acquire();
// auto stencil_cell_list = stencil[cell_j_id];
// ShrinkMatrixView B(B_pool[t], stencil_cell_list.size());
// size_t column_begin = 0;
// for (size_t i_discrete_function_variant = 0;
// i_discrete_function_variant < discrete_function_variant_list.size(); ++i_discrete_function_variant) {
// // const auto& discrete_function_variant = discrete_function_variant_list[i_discrete_function_variant];
// // std::visit(
// // [&](auto&& discrete_function) {
// // using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
// // if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
// // using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
// // const DataType& pj = discrete_function[cell_j_id];
// // for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
// // const CellId cell_i_id = stencil_cell_list[i];
// // const DataType& pi_pj = discrete_function[cell_i_id] - pj;
// // if constexpr (std::is_arithmetic_v<DataType>) {
// // B(i, column_begin) = pi_pj;
// // } else if constexpr (is_tiny_vector_v<DataType>) {
// // for (size_t kB = column_begin, k = 0; k < DataType::Dimension; ++k, ++kB) {
// // B(i, kB) = pi_pj[k];
// // }
// // } else if constexpr (is_tiny_matrix_v<DataType>) {
// // for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
// // for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
// // B(i, column_begin + k * DataType::NumberOfColumns + l) = pi_pj(k, l);
// // }
// // }
// // }
// // }
// // if constexpr (std::is_arithmetic_v<DataType>) {
// // ++column_begin;
// // } else if constexpr (is_tiny_vector_v<DataType> or is_tiny_matrix_v<DataType>) { // // column_begin += DataType::Dimension;
// // }
// // }
// // },
// // discrete_function_variant->discreteFunction());
// using DataType = double;
// const DataType& pj = my_discrete_function[cell_j_id];
// for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
// const CellId cell_i_id = stencil_cell_list[i];
// const DataType& pi_pj = my_discrete_function[cell_i_id] - pj;
// B(i, column_begin) = pi_pj;
// }
// ++column_begin;
// }
// ShrinkMatrixView A(A_pool[t], stencil_cell_list.size());
// const Rd& Xj = xj[cell_j_id];
// for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
// const CellId cell_i_id = stencil_cell_list[i];
// const Rd Xi_Xj = xj[cell_i_id] - Xj;
// for (size_t l = 0; l < MeshType::Dimension; ++l) {
// A(i, l) = Xi_Xj[l];
// }
// }
// const SmallMatrix<double>& X = X_pool[t];
// Givens::solveCollectionInPlace(A, X, B);
// column_begin = 0;
// for (size_t i_discrete_function_variant = 0;
// i_discrete_function_variant < discrete_function_variant_list.size(); ++i_discrete_function_variant) {
// auto discrete_function_variant = discrete_function_variant_list[i_discrete_function_variant];
// std::visit(
// [&](auto&& discrete_function) {
// using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
// if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
// using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
// const DataType& pj = discrete_function[cell_j_id];
// auto discrete_function_dpk = mutable_discrete_function_dpk_variant_list[i_discrete_function_variant]
// .get<DiscreteFunctionDPk<MeshType::Dimension, DataType>>();
// auto dpk_j = discrete_function_dpk.coefficients(cell_j_id);
// dpk_j[0] = pj;
// if constexpr (std::is_arithmetic_v<DataType>) {
// for (size_t i = 0; i < MeshType::Dimension; ++i) {
// auto& dpk_j_ip1 = dpk_j[i + 1];
// dpk_j_ip1 = X(i, column_begin);
// }
// ++column_begin;
// } else if constexpr (is_tiny_vector_v<DataType>) {
// for (size_t i = 0; i < MeshType::Dimension; ++i) {
// auto& dpk_j_ip1 = dpk_j[i + 1];
// for (size_t k = 0; k < DataType::Dimension; ++k) {
// dpk_j_ip1[k] = X(i, column_begin + k);
// }
// }
// column_begin += DataType::Dimension;
// } else if constexpr (is_tiny_matrix_v<DataType>) {
// for (size_t i = 0; i < MeshType::Dimension; ++i) {
// auto& dpk_j_ip1 = dpk_j[i + 1];
// for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
// for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
// dpk_j_ip1(k, l) = X(i, column_begin + k * DataType::NumberOfColumns + l);
// } // }
// }
// column_begin += DataType::Dimension;
// }
// } else {
// throw NotImplementedError("discrete function type");
// }
// },
// discrete_function_variant->discreteFunction());
// }
// tokens.release(t);
// }
// });
std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> discrete_function_dpk_variant_list;
// for (auto discrete_function_dpk_variant_p : mutable_discrete_function_dpk_variant_list) {
// discrete_function_dpk_variant_list.push_back(
// std::make_shared<const DiscreteFunctionDPkVariant>(discrete_function_dpk_variant_p));
// }
return discrete_function_dpk_variant_list;
}
std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
PolynomialReconstruction::build(
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
{
if (not hasSameMesh(discrete_function_variant_list)) {
throw NormalError("cannot reconstruct functions living of different mesh simultaneously");
}
auto mesh_v = getCommonMesh(discrete_function_variant_list);
return std::visit([&](auto&& p_mesh) { return this->_build(p_mesh, discrete_function_variant_list); },
mesh_v->variant());
}