diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp index 0be440945d5b9f3763837192e3e4f1ff0eaf3f5d..76dfefd4a791db45b2f5cb565df4a241b14bfecd 100644 --- a/src/language/modules/SchemeModule.cpp +++ b/src/language/modules/SchemeModule.cpp @@ -425,12 +425,12 @@ SchemeModule::SchemeModule() std::function( [](const std::shared_ptr<const DiscreteFunctionVariant>& rho, - const std::shared_ptr<const DiscreteFunctionVariant>& U, - const std::shared_ptr<const DiscreteFunctionVariant>& E, + const std::shared_ptr<const DiscreteFunctionVariant>& rho_U, + const std::shared_ptr<const DiscreteFunctionVariant>& rho_E, const double gamma, 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); + return gks(rho, rho_U, rho_E, gamma, dt); } )); diff --git a/src/scheme/GKS.cpp b/src/scheme/GKS.cpp index 09fa5f4ff893cd6a37847b2fb07f7c246e19cdf1..c3f296b73d938b8e520e83892d4a91b143a8e20b 100644 --- a/src/scheme/GKS.cpp +++ b/src/scheme/GKS.cpp @@ -15,43 +15,125 @@ class GKS 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, + std::shared_ptr<const DiscreteFunctionVariant> rho_U_v, + std::shared_ptr<const DiscreteFunctionVariant> rho_E_v, + const double gamma, 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>>(); + const double pi = acos(-1); + const int delta = (3 - gamma) / (gamma - 1); - DiscreteFunctionP0<double> rho = copy(rho_n); - DiscreteFunctionP0<Rd> U = copy(U_n); - DiscreteFunctionP0<double> E = copy(E_n); + DiscreteFunctionP0<const double> rho_n = rho_v->get<DiscreteFunctionP0<const double>>(); + DiscreteFunctionP0<const Rd> rho_U_n = rho_U_v->get<DiscreteFunctionP0<const Rd>>(); + DiscreteFunctionP0<const double> rho_E_n = rho_E_v->get<DiscreteFunctionP0<const double>>(); + + DiscreteFunctionP0<double> rho = copy(rho_n); + DiscreteFunctionP0<Rd> rho_U = copy(rho_U_n); + DiscreteFunctionP0<double> rho_E = copy(rho_E_n); + + CellValue<double> lambda{p_mesh->connectivity()}; + CellValue<Rd> U{p_mesh->connectivity()}; auto& mesh_data = MeshDataManager::instance().getMeshData(mesh); auto Vj = mesh_data.Vj(); auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix(); + // const NodeId first_node_id = NodeId{0}; + // const NodeId last_node_id = NodeId{static_cast<typename CellId::base_type>(mesh.numberOfNodes() - 1)}; NodeValue<double> rho_flux(mesh.connectivity()); + NodeValue<Rd> rho_U_flux(mesh.connectivity()); + NodeValue<double> rho_E_flux(mesh.connectivity()); + + NodeValue<double> rho_node(mesh.connectivity()); + NodeValue<Rd> rho_U_node(mesh.connectivity()); + NodeValue<double> rho_E_node(mesh.connectivity()); + + for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) { + U[cell_id][0] = rho_U_n[cell_id][0] / rho_n[cell_id]; + double U_2 = U[cell_id][0] * U[cell_id][0]; + double rho_U_2 = rho_U_n[cell_id][0] * U[cell_id][0]; + lambda[cell_id] = (1. + delta) * (rho_n[cell_id]) / (4 * rho_E_n[cell_id] - 2 * rho_U_2); + + double rho_cell_left = rho_n[cell_id] * (1 + std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0])); + rho_cell_left /= 2.; + + Rd rho_U_cell_left; + rho_U_cell_left[0] = rho_U_n[cell_id][0] * (1. + std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0])) + + rho_n[cell_id] * std::exp(-lambda[cell_id] * U_2) / std::sqrt(pi * lambda[cell_id]); + rho_U_cell_left[0] /= 2.; + + double rho_E_cell_left = + rho_E_n[cell_id] * (1. + std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0])) + + 0.5 * rho_U_n[cell_id][0] * std::exp(-lambda[cell_id] * U_2) / std::sqrt(pi * lambda[cell_id]); + rho_E_cell_left /= 2.; + + auto node_list = cell_to_node_matrix[cell_id]; + + rho_node[node_list[1]] = rho_cell_left; + rho_U_node[node_list[1]] = rho_U_cell_left; + rho_E_node[node_list[1]] = rho_E_cell_left; + } + 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]; + double U_2 = U[cell_id][0] * U[cell_id][0]; + + double rho_cell_right = rho_n[cell_id] * (1 - std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0])); + rho_cell_right /= 2.; + + Rd rho_U_cell_right; + rho_U_cell_right[0] = rho_U_n[cell_id][0] * (1. - std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0])) - + rho_n[cell_id] * std::exp(-lambda[cell_id] * U_2) / std::sqrt(pi * lambda[cell_id]); + rho_U_cell_right[0] /= 2.; + + double rho_E_cell_right = + rho_E_n[cell_id] * (1. - std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0])) - + 0.5 * rho_U_n[cell_id][0] * std::exp(-lambda[cell_id] * U_2) / std::sqrt(pi * lambda[cell_id]); + rho_E_cell_right /= 2.; + + auto node_list = cell_to_node_matrix[cell_id]; + + rho_node[node_list[0]] += rho_cell_right; + rho_U_node[node_list[0]] += rho_U_cell_right; + rho_E_node[node_list[0]] += rho_E_cell_right; } - 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]; + auto node_list = cell_to_node_matrix[cell_id]; + + rho_flux[node_list[1]] = rho_U_node[node_list[1]][0]; + rho_U_flux[node_list[1]][0] = + delta / (1. + delta) * rho_U_node[node_list[1]][0] * rho_U_node[node_list[1]][0] / rho_node[node_list[1]] + + 2. / (1. + delta) * rho_E_node[node_list[1]]; + rho_E_flux[node_list[1]] = + rho_U_node[node_list[1]][0] / rho_node[node_list[1]] * + ((3. + delta) / (1. + delta) * rho_E_node[node_list[1]] - + 1. / (1. + delta) * rho_U_node[node_list[1]][0] * rho_U_node[node_list[1]][0] / rho_node[node_list[1]]); + } + + // Peut être un problème ici avec les bords mais je sais pas encore + // rho_flux[first_node_id] = rho_flux[last_node_id]; + // rho_U_flux[first_node_id] = rho_U_flux[last_node_id]; + // rho_E_flux[first_node_id] = rho_E_flux[last_node_id]; + + for (CellId cell_id = 1; cell_id < mesh.numberOfCells() - 1; ++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]]); + const Rd rho_U_flux_sum = (rho_U_flux[node_list[1]] - rho_U_flux[node_list[0]]); + const double rho_E_flux_sum = (rho_E_flux[node_list[1]] - rho_E_flux[node_list[0]]); - rho[cell_id] -= dt * 0.5 * rho_flux_sum; + rho[cell_id] -= dt * rho_flux_sum / Vj[cell_id]; + rho_U[cell_id][0] -= dt * rho_U_flux_sum[0] / Vj[cell_id]; + rho_E[cell_id] -= dt * rho_E_flux_sum / Vj[cell_id]; } - return std::make_tuple(std::make_shared<DiscreteFunctionVariant>(rho), std::make_shared<DiscreteFunctionVariant>(U), - std::make_shared<DiscreteFunctionVariant>(E)); + return std::make_tuple(std::make_shared<DiscreteFunctionVariant>(rho), + std::make_shared<DiscreteFunctionVariant>(rho_U), + std::make_shared<DiscreteFunctionVariant>(rho_E)); } GKS() = default; @@ -61,16 +143,17 @@ 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, + std::shared_ptr<const DiscreteFunctionVariant> rho_U_v, + std::shared_ptr<const DiscreteFunctionVariant> rho_E_v, + const double gamma, const double dt) { - std::shared_ptr mesh_v = getCommonMesh({rho_v, U_v, E_v}); + std::shared_ptr mesh_v = getCommonMesh({rho_v, rho_U_v, rho_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)) { + if (not checkDiscretizationType({rho_v, rho_U_v, rho_E_v}, DiscreteFunctionType::P0)) { throw NormalError("GKS solver expects P0 functions"); } @@ -82,7 +165,7 @@ gks(std::shared_ptr<const DiscreteFunctionVariant> rho_v, 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); + return gks.solve(p_mesh, rho_v, rho_U_v, rho_E_v, gamma, dt); } else { throw NormalError("dimension not treated"); diff --git a/src/scheme/GKS.hpp b/src/scheme/GKS.hpp index cd3643003a9e7b7be8391a5d41e3d0607c34bce5..e29ad2aaa49414c3db84e2c14e78f22cd0c66679 100644 --- a/src/scheme/GKS.hpp +++ b/src/scheme/GKS.hpp @@ -11,6 +11,7 @@ std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, // rho gks(std::shared_ptr<const DiscreteFunctionVariant> rho, std::shared_ptr<const DiscreteFunctionVariant> U, std::shared_ptr<const DiscreteFunctionVariant> E, + const double gamma, const double dt); #endif // GKS_HPP