From 04e4f9fcda3838770839b10ff387a72a580b5216 Mon Sep 17 00:00:00 2001 From: LABOURASSE Emmanuel <labourassee@gmail.com> Date: Tue, 15 Sep 2020 17:21:12 +0200 Subject: [PATCH] add the files to use Discontinuous Galerkin for 1D ODE --- .../algorithms/ODEDiscontinuousGalerkin1D.cpp | 86 +++++++ .../algorithms/ODEDiscontinuousGalerkin1D.hpp | 16 ++ src/language/utils/InterpolateArray.hpp | 44 ++++ src/scheme/DiscontinuousGalerkin1DTools.hpp | 77 +++++++ src/scheme/ODEGD1D.hpp | 214 ++++++++++++++++++ tests/test_DiscontinuousGalerkin1D.cpp | 56 +++++ 6 files changed, 493 insertions(+) create mode 100644 src/language/algorithms/ODEDiscontinuousGalerkin1D.cpp create mode 100644 src/language/algorithms/ODEDiscontinuousGalerkin1D.hpp create mode 100644 src/language/utils/InterpolateArray.hpp create mode 100644 src/scheme/DiscontinuousGalerkin1DTools.hpp create mode 100644 src/scheme/ODEGD1D.hpp create mode 100644 tests/test_DiscontinuousGalerkin1D.cpp diff --git a/src/language/algorithms/ODEDiscontinuousGalerkin1D.cpp b/src/language/algorithms/ODEDiscontinuousGalerkin1D.cpp new file mode 100644 index 000000000..4bd10ddc0 --- /dev/null +++ b/src/language/algorithms/ODEDiscontinuousGalerkin1D.cpp @@ -0,0 +1,86 @@ +#include <language/algorithms/ODEDiscontinuousGalerkin1D.hpp> +#include <language/utils/InterpolateArray.hpp> +#include <scheme/ODEGD1D.hpp> + +#include <fstream> + +template <size_t Dimension, size_t Degree> +ODEDiscontinuousGalerkin1D<Dimension, Degree>::ODEDiscontinuousGalerkin1D(const BasisType& basis_type, + std::shared_ptr<const IMesh> i_mesh, + const FunctionSymbolId& uex_id) +{ + using ConnectivityType = Connectivity<Dimension>; + using MeshType = Mesh<ConnectivityType>; + + std::shared_ptr<const MeshType> mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh); + static_assert(Dimension == 1, "Dimension have to be 1 for DiscontinuousGalerkin1D"); + using R1 = TinyVector<Dimension>; + ODEDG1D<Degree> odedg1d(basis_type, mesh); + const NodeValue<const R1> xr = mesh->xr(); + const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix(); + Array<double> my_array; + for (CellId j = 0; j < mesh->numberOfCells(); ++j) { + const auto& cell_nodes = cell_to_node_matrix[j]; + const NodeId r1 = cell_nodes[0]; + const NodeId r2 = cell_nodes[1]; + const TinyVector<Dimension> xr1 = xr[r1]; + const TinyVector<Dimension> xr2 = xr[r2]; + double x13 = (2 * xr1[0] + xr2[0]) / 3; + double x23 = (xr1[0] + 2 * xr2[0]) / 3; + my_array[j * 4] = xr1[0]; + my_array[j * 4 + 1] = x13; + my_array[j * 4 + 2] = x23; + my_array[j * 4 + 3] = xr2[0]; + } + + Array<double> ej = InterpolateArray<double(double)>::interpolate(uex_id, my_array); + + // InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(uex_id, mesh_data.xj()); + + CellValue<Polynomial<Degree>> U(mesh->connectivity()); + odedg1d.globalSolve(U); + PolynomialBasis<Degree> B; + TinyVector<Degree + 1> zeros; + for (size_t i = 0; i <= Degree; i++) { + zeros[i] = i; + } + B.build(basis_type, 0.5, zeros); + std::string filename = "result-basis-"; + filename += B.displayType(); + filename += "-degree-"; + filename += std::to_string(Degree); + filename += "-dof-"; + filename += std::to_string(mesh->numberOfCells()); + std::string filename2 = filename + "-polynomials"; + std::ofstream sout(filename.c_str()); + std::ofstream sout2(filename2.c_str()); + for (CellId j = 0; j < mesh->numberOfCells(); ++j) { + const auto& cell_nodes = cell_to_node_matrix[j]; + const NodeId r1 = cell_nodes[0]; + const NodeId r2 = cell_nodes[1]; + const TinyVector<Dimension> xr1 = xr[r1]; + const TinyVector<Dimension> xr2 = xr[r2]; + const double deltax = (xr2[0] - xr1[0]) / 100; + for (size_t i = 0; i <= 100; ++i) { + const double x = xr1[0] + deltax * i; + sout2 << x << " " << U[j](x) << '\n'; + } + double x13 = (2 * xr1[0] + xr2[0]) / 3; + double x23 = (xr1[0] + 2 * xr2[0]) / 3; + + double uex = std::exp(xr1[0]) + 3 * (std::exp(x13) + std::exp(x23)) + std::exp(xr2[0]); + uex /= 8; + double ucalc = integrate(U[j], xr1[0], xr2[0]) / (xr2[0] - xr1[0]); + sout << (0.5 * (xr1[0] + xr2[0])) << " " << uex << " " << ucalc << " " << std::abs(ucalc - uex) << '\n'; + } +} + +template ODEDiscontinuousGalerkin1D<1, 0>::ODEDiscontinuousGalerkin1D(const BasisType& basis_type, + std::shared_ptr<const IMesh>, + const FunctionSymbolId&); +template ODEDiscontinuousGalerkin1D<1, 1>::ODEDiscontinuousGalerkin1D(const BasisType& basis_type, + std::shared_ptr<const IMesh>, + const FunctionSymbolId&); +template ODEDiscontinuousGalerkin1D<1, 2>::ODEDiscontinuousGalerkin1D(const BasisType& basis_type, + std::shared_ptr<const IMesh>, + const FunctionSymbolId&); diff --git a/src/language/algorithms/ODEDiscontinuousGalerkin1D.hpp b/src/language/algorithms/ODEDiscontinuousGalerkin1D.hpp new file mode 100644 index 000000000..f7425d8d9 --- /dev/null +++ b/src/language/algorithms/ODEDiscontinuousGalerkin1D.hpp @@ -0,0 +1,16 @@ +#ifndef ORDINARY_DIFFERENTIAL_EQUATION_DISCONTINUOUS_GALERKIN_1D_HPP +#define ORDINARY_DIFFERENTIAL_EQUATION_DISCONTINUOUS_GALERKIN_1D_HPP + +#include <analysis/PolynomialBasis.hpp> +#include <language/utils/FunctionSymbolId.hpp> +#include <mesh/IMesh.hpp> + +template <size_t Dimension, size_t Degree> +struct ODEDiscontinuousGalerkin1D +{ + ODEDiscontinuousGalerkin1D(const BasisType& basis_type, + std::shared_ptr<const IMesh> i_mesh, + const FunctionSymbolId& uex_id); +}; + +#endif // ORDINARY_DIFFERENTIAL_EQUATION_DISCONTINUOUS_GALERKIN_1D_HPP diff --git a/src/language/utils/InterpolateArray.hpp b/src/language/utils/InterpolateArray.hpp new file mode 100644 index 000000000..be5c114f5 --- /dev/null +++ b/src/language/utils/InterpolateArray.hpp @@ -0,0 +1,44 @@ +#ifndef INTERPOLATE_ARRAY_HPP +#define INTERPOLATE_ARRAY_HPP + +#include <language/utils/PugsFunctionAdapter.hpp> + +template <typename T> +class InterpolateArray; +template <typename OutputType, typename InputType> +class InterpolateArray<OutputType(InputType)> : public PugsFunctionAdapter<OutputType(InputType)> +{ + static constexpr size_t Dimension = OutputType::Dimension; + using Adapter = PugsFunctionAdapter<OutputType(InputType)>; + + public: + static inline Array<OutputType> + interpolate(const FunctionSymbolId& function_symbol_id, const Array<const InputType>& position) + { + auto& expression = Adapter::getFunctionExpression(function_symbol_id); + auto convert_result = Adapter::getResultConverter(expression.m_data_type); + + Array<ExecutionPolicy> context_list = Adapter::getContextList(expression); + + using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space; + Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens; + + Array<OutputType> value(position.size()); + + parallel_for(position.size(), [=, &expression, &tokens](size_t i) { + const int32_t t = tokens.acquire(); + + auto& execution_policy = context_list[t]; + + Adapter::convertArgs(execution_policy.currentContext(), position[i]); + auto result = expression.execute(execution_policy); + value[i] = convert_result(std::move(result)); + + tokens.release(t); + }); + + return value; + } +}; + +#endif /* INTERPOLATE_ARRAY_HPP */ diff --git a/src/scheme/DiscontinuousGalerkin1DTools.hpp b/src/scheme/DiscontinuousGalerkin1DTools.hpp new file mode 100644 index 000000000..d28a03370 --- /dev/null +++ b/src/scheme/DiscontinuousGalerkin1DTools.hpp @@ -0,0 +1,77 @@ +#ifndef DISCONTINUOUS_GALERKIN_1D_TOOLS +#define DISCONTINUOUS_GALERKIN_1D_TOOLS + +#include <rang.hpp> + +#include <utils/ArrayUtils.hpp> + +#include <utils/PugsAssert.hpp> + +#include <mesh/Mesh.hpp> +#include <mesh/MeshData.hpp> +#include <mesh/MeshDataManager.hpp> +#include <mesh/MeshNodeBoundary.hpp> + +#include <algebra/TinyMatrix.hpp> +#include <algebra/TinyVector.hpp> + +#include <analysis/Polynomial.hpp> +#include <analysis/PolynomialBasis.hpp> + +#include <mesh/ItemValueUtils.hpp> +#include <mesh/SubItemValuePerItem.hpp> + +#include <utils/Exceptions.hpp> +#include <utils/Messenger.hpp> + +#include <iostream> + +template <size_t Degree> +class DiscontinuousGalerkin1DTools +{ + constexpr static size_t Dimension = 1; + + using Rd = TinyVector<Degree>; + using Rdd = TinyMatrix<Degree>; + + using R1 = TinyVector<Dimension>; + + public: + int + inverse() const + { + return 0; + } + + void + elementarySolve(Polynomial<Degree>& U, + const TinyVector<Degree + 1> moments, + const PolynomialBasis<Degree> Basis, + const double& xinf, + const double& xsup) const + { + TinyMatrix<Degree + 1> M(zero); + for (size_t i = 0; i <= Degree; ++i) { + for (size_t j = 0; j <= Degree; ++j) { + Polynomial<2 * Degree> P = Basis.elements()[i] * Basis.elements()[j]; + M(i, j) = integrate(P, xinf, xsup); + } + } + TinyMatrix inv_M = ::inverse(M); + U.coefficients() = inv_M * moments; + } + + // void + // integrate() const + // {} + + public: + // TODO add a flux manager to constructor + // TODO add a basis to constructor + DiscontinuousGalerkin1DTools() + { + ; + } +}; + +#endif // DISCONTINUOUS_GALERKIN_1D_TOOLS diff --git a/src/scheme/ODEGD1D.hpp b/src/scheme/ODEGD1D.hpp new file mode 100644 index 000000000..c2d89c22b --- /dev/null +++ b/src/scheme/ODEGD1D.hpp @@ -0,0 +1,214 @@ +#ifndef ODE_DG_1D +#define ODE_DG_1D + +#include <rang.hpp> + +#include <utils/ArrayUtils.hpp> + +#include <utils/PugsAssert.hpp> + +#include <mesh/Mesh.hpp> +#include <mesh/MeshData.hpp> +#include <mesh/MeshDataManager.hpp> +#include <mesh/MeshNodeBoundary.hpp> + +#include <algebra/TinyMatrix.hpp> +#include <algebra/TinyVector.hpp> + +#include <analysis/Polynomial.hpp> +#include <analysis/PolynomialBasis.hpp> + +#include <mesh/ItemValueUtils.hpp> +#include <mesh/SubItemValuePerItem.hpp> + +#include <scheme/DiscontinuousGalerkin1DTools.hpp> +#include <utils/Exceptions.hpp> +#include <utils/Messenger.hpp> + +#include <iostream> +enum class FluxType +{ + centered, + stable +}; + +template <size_t Degree> +class ODEDG1D +{ + using MeshType = Mesh<Connectivity1D>; + constexpr static size_t Dimension = MeshType::Dimension; + static_assert(Dimension == 1, "Galerkin 1D only works in dimension 1"); + using MeshDataType = MeshData<Dimension>; + // using UnknownsType = FiniteVolumesEulerUnknowns<MeshType>; + + const BasisType& m_basis_type; + std::shared_ptr<const MeshType> m_mesh; + const typename MeshType::Connectivity& m_connectivity; + + // class DirichletBoundaryCondition; + + // using BoundaryCondition = std::variant<DirichletBoundaryCondition>; + + // using BoundaryConditionList = std::vector<BoundaryCondition>; + // DiscontinuousGalerkin1DTools<Degree>& m_dgtools; + using Rd = TinyVector<Degree>; + using Rdd = TinyMatrix<Degree>; + + PolynomialBasis<Degree> m_basis; + using R1 = TinyVector<Dimension>; + // const BoundaryConditionList m_boundary_condition_list; + // NodeValue<R1> m_flux; + TinyVector<Degree + 1> + _buildZeros(double xmin, double xmax) + { + TinyVector<Degree + 1> zeros(zero); + zeros[0] = xmin; + double dx = (xmax - xmin) / Degree; + for (size_t i = 1; i <= Degree; ++i) { + zeros[i] = xmin + i * dx; + } + return zeros; + } + + public: + // void + // _applyBC() + // { + // for (const auto& boundary_condition : m_boundary_condition_list) { + // std::visit( + // [&](auto&& bc) { + // using T = std::decay_t<decltype(bc)>; + // if constexpr (std::is_same_v<DirichletBoundaryCondition, T>) { + // const auto& node_list = bc.nodeList(); + // const auto& value_list = bc.valueList(); + + // parallel_for( + // node_list.size(), PUGS_LAMBDA(size_t i_node) { + // NodeId node_id = node_list[i_node]; + // const auto& value = value_list[i_node]; + + // m_flux[node_id] = value; + // }); + // } + // }, + // boundary_condition); + // } + // } + + // NodeValue<R1> + // computeFluxes(FluxType flux_type, const PolynomialBasis<Degree> Basis, Polynomial<Degree>& U) + // { + // switch (flux_type) { + // // case FluxType::centered: { + // // return computeCenteredFlux(Basis, U); + // // break; + // // } + // case FluxType::stable: { + // return computeStableFlux(Basis, U); + // break; + // } + // default: + // throw UnexpectedError("unknown flux type"); + // } + // } + + double + computeStableFluxes(const CellValue<const Polynomial<Degree>>& U, NodeId r) + { + double flux; + const NodeValue<const R1> xr = m_mesh->xr(); + const auto& node_to_cell_matrix = m_mesh->connectivity().nodeToCellMatrix(); + + const auto& node_cells = node_to_cell_matrix[r]; + if (node_cells.size() == 1) { + flux = 1.; + } else { + const CellId j1 = node_cells[0]; + const double xr0 = xr[r][0]; + flux = U[j1](xr0); + std::cout << "U[" << j1 << "]=" << U[j1] << "\n"; + } + return flux; + } + + void + globalSolve(CellValue<Polynomial<Degree>>& U) + { + const NodeValue<const R1> xr = m_mesh->xr(); + const auto& cell_to_node_matrix = m_mesh->connectivity().cellToNodeMatrix(); + CellValue<TinyVector<Degree + 1>> moments(m_connectivity); + for (CellId j = 0; j < m_mesh->numberOfCells(); j++) { + const auto& cell_nodes = cell_to_node_matrix[j]; + const NodeId r1 = cell_nodes[0]; + const NodeId r2 = cell_nodes[1]; + const TinyVector<Dimension> xr1 = xr[r1]; + const TinyVector<Dimension> xr2 = xr[r2]; + const TinyVector<Degree + 1> zeros = _buildZeros(xr1[0], xr2[0]); + std::cout << "x1 " << xr1[0] << " x2 " << xr2[0] << " zeros " << zeros << "\n"; + m_basis.build(m_basis_type, 0.5 * (xr1[0] + xr2[0]), zeros); + std::cout << " B[" << j << "]=" << m_basis.elements()[0] << "\n"; + const double fluxg = computeStableFluxes(U, r1); + std::cout << j << " flux " << fluxg << "\n"; + moments[j] = buildMoments(m_basis, fluxg, xr1[0]); + std::cout << " moments[" << j << "]=" << moments[j] << "\n"; + + U[j] = elementarySolve(moments[j], m_basis, xr1[0], xr2[0]); + std::cout << " U[" << j << "]=" << U[j] << "\n"; + } + } + + // void + // integrate() const + // {} + PUGS_INLINE constexpr TinyVector<Degree + 1> + buildMoments(PolynomialBasis<Degree> basis, double fluxg, double xr1) + { + TinyVector<Degree + 1> moments(zero); + for (size_t i = 0; i <= Degree; i++) { + // moments[i] = (basis.elements()[i](xr2) - basis.elements()[i](xr1)) * fluxg; + moments[i] = -basis.elements()[i](xr1) * fluxg; + } + return moments; + } + + PUGS_INLINE constexpr Polynomial<Degree> + elementarySolve(const TinyVector<Degree + 1> moments, + const PolynomialBasis<Degree> Basis, + const double& xinf, + const double& xsup) const + { + Polynomial<Degree> U; + U.coefficients() = zero; + TinyMatrix<Degree + 1> M(zero); + for (size_t i = 0; i <= Degree; ++i) { + for (size_t j = 0; j <= Degree; ++j) { + Polynomial<2 * Degree> P = Basis.elements()[j] * Basis.elements()[i]; + P += derivative(Basis.elements()[i]) * Basis.elements()[j]; + double fluxd = Basis.elements()[i](xsup) * Basis.elements()[j](xsup); + M(i, j) = integrate(P, xinf, xsup) - fluxd; + } + } + std::cout << " M " << M << "\n"; + TinyMatrix inv_M = ::inverse(M); + TinyVector<Degree + 1> coefficients = inv_M * moments; + for (size_t i = 0; i <= Degree; ++i) { + U += coefficients[i] * Basis.elements()[i]; + } + std::cout << "U(" << xsup << ")=" << U(xsup) << "\n"; + return U; + } + + public: + // TODO add a flux manager to constructor + // TODO add a basis to constructor + ODEDG1D(const BasisType& basis_type, std::shared_ptr<const IMesh> i_mesh) //, const BoundaryConditionList& bc_list) + : m_basis_type(basis_type), + m_mesh(std::dynamic_pointer_cast<const MeshType>(i_mesh)), + m_connectivity(m_mesh->connectivity()) //, + // m_boundary_condition_list(bc_list) + { + ; + } +}; + +#endif // ODE_DG_1D diff --git a/tests/test_DiscontinuousGalerkin1D.cpp b/tests/test_DiscontinuousGalerkin1D.cpp new file mode 100644 index 000000000..cb4109d23 --- /dev/null +++ b/tests/test_DiscontinuousGalerkin1D.cpp @@ -0,0 +1,56 @@ +#include <catch2/catch.hpp> + +#include <Kokkos_Core.hpp> + +#include <utils/PugsAssert.hpp> +#include <utils/Types.hpp> + +#include <algebra/TinyMatrix.hpp> +#include <analysis/Polynomial.hpp> +#include <analysis/PolynomialBasis.hpp> +#include <mesh/CartesianMeshBuilder.hpp> +#include <mesh/DiamondDualConnectivityManager.hpp> +#include <mesh/DiamondDualMeshManager.hpp> +#include <scheme/DiscontinuousGalerkin1DTools.hpp> +#include <scheme/ODEGD1D.hpp> + +// Instantiate to ensure full coverage is performed +template class Polynomial<0>; +template class Polynomial<1>; +template class Polynomial<2>; +template class Polynomial<3>; +template class Polynomial<4>; +template class Polynomial<5>; + +// clazy:excludeall=non-pod-global-static + +TEST_CASE("Discontinuous-Galerkin-1D", "[discontinuous-galerkin-1d]") +{ + SECTION("Elementary-tests") + { + TinyVector<1> a{0}; + TinyVector<1> b{1}; + TinyVector<1, uint64_t> nbcells{10}; + // CartesianMeshBuilder mesh_builder = CartesianMeshBuilder(a, b, nbcells); + // auto mesh = mesh_builder.mesh(); + // DiscontinuousGalerkin1D<2>{mesh}; + REQUIRE_NOTHROW(DiscontinuousGalerkin1DTools<2>()); + DiscontinuousGalerkin1DTools dg = DiscontinuousGalerkin1DTools<1>(); + Polynomial<1> U; + TinyVector<2> F({3, 2}); + PolynomialBasis<1> Basis; + Basis.build(BasisType::canonical); + dg.elementarySolve(U, F, Basis, -1, 1); + REQUIRE(U == Polynomial<1>{{3. / 2, 3}}); + Polynomial<2> U2; + TinyVector<3> F2({1, 1, 1}); + PolynomialBasis<2> Basis2; + Basis2.build(BasisType::canonical); + DiscontinuousGalerkin1DTools dg2 = DiscontinuousGalerkin1DTools<2>(); + dg2.elementarySolve(U2, F2, Basis2, -1, 1); + REQUIRE(integrate(U2, -1, 1) == Approx(1).epsilon(1E-14)); + REQUIRE(integrate(U2 * Polynomial<2>{{0, 0, 1}}, -1, 1) == 1.); + dg2.elementarySolve(U2, F2, Basis2, 0, 1); + REQUIRE(integrate(U2 * Polynomial<2>{{0, 0, 1}}, 0, 1) == Approx(1).epsilon(1E-14)); + } +} -- GitLab