#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); if constexpr (Degree == 0) { zeros[0] = 0.5 * (xmax - xmin); } else { 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