diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp index 7b404fcf670c5217721efba890728ce97e3e91a4..0be440945d5b9f3763837192e3e4f1ff0eaf3f5d 100644 --- a/src/language/modules/SchemeModule.cpp +++ b/src/language/modules/SchemeModule.cpp @@ -35,6 +35,7 @@ #include <scheme/FluxingAdvectionSolver.hpp> #include <scheme/FourierBoundaryConditionDescriptor.hpp> #include <scheme/FreeBoundaryConditionDescriptor.hpp> +#include <scheme/GKS.hpp> #include <scheme/HyperelasticSolver.hpp> #include <scheme/IBoundaryConditionDescriptor.hpp> #include <scheme/IDiscreteFunctionDescriptor.hpp> @@ -420,6 +421,20 @@ SchemeModule::SchemeModule() )); + this->_addBuiltinFunction("gks", + std::function( + + [](const std::shared_ptr<const DiscreteFunctionVariant>& rho, + const std::shared_ptr<const DiscreteFunctionVariant>& U, + const std::shared_ptr<const DiscreteFunctionVariant>& E, + const double& dt) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, + std::shared_ptr<const DiscreteFunctionVariant>, + std::shared_ptr<const DiscreteFunctionVariant>> { + return gks(rho, U, E, dt); + } + + )); + this->_addBuiltinFunction("glace_solver", std::function( diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt index dadbe73bb24ad84e491248b89016baf85f0dbc31..f3aa4dc1b0704ebd319d767705582f74e87d9464 100644 --- a/src/scheme/CMakeLists.txt +++ b/src/scheme/CMakeLists.txt @@ -10,7 +10,9 @@ add_library( DiscreteFunctionVectorIntegrator.cpp DiscreteFunctionVectorInterpoler.cpp FluxingAdvectionSolver.cpp -) + HyperelasticSolver.cpp + GKS.cpp + ) target_link_libraries( PugsScheme diff --git a/src/scheme/GKS.cpp b/src/scheme/GKS.cpp new file mode 100644 index 0000000000000000000000000000000000000000..09fa5f4ff893cd6a37847b2fb07f7c246e19cdf1 --- /dev/null +++ b/src/scheme/GKS.cpp @@ -0,0 +1,96 @@ +#include <scheme/GKS.hpp> + +#include <mesh/Mesh.hpp> +#include <mesh/MeshData.hpp> +#include <mesh/MeshDataManager.hpp> +#include <mesh/MeshTraits.hpp> +#include <scheme/DiscreteFunctionUtils.hpp> + +template <MeshConcept MeshType> +class GKS +{ + public: + std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, + std::shared_ptr<const DiscreteFunctionVariant>, + std::shared_ptr<const DiscreteFunctionVariant>> + solve(std::shared_ptr<const MeshType> p_mesh, + std::shared_ptr<const DiscreteFunctionVariant> rho_v, + std::shared_ptr<const DiscreteFunctionVariant> U_v, + std::shared_ptr<const DiscreteFunctionVariant> E_v, + double dt) + { + using Rd = TinyVector<MeshType::Dimension>; + + const MeshType& mesh = *p_mesh; + + DiscreteFunctionP0<const double> rho_n = rho_v->get<DiscreteFunctionP0<const double>>(); + DiscreteFunctionP0<const Rd> U_n = U_v->get<DiscreteFunctionP0<const Rd>>(); + DiscreteFunctionP0<const double> E_n = E_v->get<DiscreteFunctionP0<const double>>(); + + DiscreteFunctionP0<double> rho = copy(rho_n); + DiscreteFunctionP0<Rd> U = copy(U_n); + DiscreteFunctionP0<double> E = copy(E_n); + + auto& mesh_data = MeshDataManager::instance().getMeshData(mesh); + auto Vj = mesh_data.Vj(); + + auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix(); + + NodeValue<double> rho_flux(mesh.connectivity()); + for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) { + auto node_list = cell_to_node_matrix[cell_id]; + rho_flux[node_list[1]] = rho_n[cell_id]; + } + rho_flux[NodeId{0}] = rho_flux[NodeId{static_cast<typename CellId::base_type>(mesh.numberOfNodes() - 1)}]; + + for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) { + auto node_list = cell_to_node_matrix[cell_id]; + const double rho_flux_sum = (rho_flux[node_list[1]] - rho_flux[node_list[0]]) / Vj[cell_id]; + + rho[cell_id] -= dt * 0.5 * rho_flux_sum; + } + + return std::make_tuple(std::make_shared<DiscreteFunctionVariant>(rho), std::make_shared<DiscreteFunctionVariant>(U), + std::make_shared<DiscreteFunctionVariant>(E)); + } + + GKS() = default; +}; + +std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, // rho + std::shared_ptr<const DiscreteFunctionVariant>, // U + std::shared_ptr<const DiscreteFunctionVariant>> // E +gks(std::shared_ptr<const DiscreteFunctionVariant> rho_v, + std::shared_ptr<const DiscreteFunctionVariant> U_v, + std::shared_ptr<const DiscreteFunctionVariant> E_v, + const double dt) +{ + std::shared_ptr mesh_v = getCommonMesh({rho_v, U_v, E_v}); + if (not mesh_v) { + throw NormalError("discrete functions are not defined on the same mesh"); + } + + if (not checkDiscretizationType({rho_v, U_v, E_v}, DiscreteFunctionType::P0)) { + throw NormalError("GKS solver expects P0 functions"); + } + + return std::visit( + [&](auto&& p_mesh) + -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>, + std::shared_ptr<const DiscreteFunctionVariant>> { + using MeshType = std::decay_t<decltype(*p_mesh)>; + if constexpr (is_polygonal_mesh_v<MeshType>) { + if constexpr (MeshType::Dimension == 1) { + GKS<MeshType> gks; + return gks.solve(p_mesh, rho_v, U_v, E_v, dt); + + } else { + throw NormalError("dimension not treated"); + } + + } else { + throw NormalError("unexpected mesh type"); + } + }, + mesh_v->variant()); +} diff --git a/src/scheme/GKS.hpp b/src/scheme/GKS.hpp new file mode 100644 index 0000000000000000000000000000000000000000..cd3643003a9e7b7be8391a5d41e3d0607c34bce5 --- /dev/null +++ b/src/scheme/GKS.hpp @@ -0,0 +1,16 @@ +#ifndef GKS_HPP +#define GKS_HPP + +#include <scheme/DiscreteFunctionVariant.hpp> + +#include <tuple> + +std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, // rho + std::shared_ptr<const DiscreteFunctionVariant>, // U + std::shared_ptr<const DiscreteFunctionVariant>> // E +gks(std::shared_ptr<const DiscreteFunctionVariant> rho, + std::shared_ptr<const DiscreteFunctionVariant> U, + std::shared_ptr<const DiscreteFunctionVariant> E, + const double dt); + +#endif // GKS_HPP