From d77b1239d7709d61ed6cf0bf0b6fdd5cf49ebfa9 Mon Sep 17 00:00:00 2001
From: Clovis <clovis.schoeck@etudiant.univ-rennes.fr>
Date: Tue, 25 Jun 2024 11:05:46 +0200
Subject: [PATCH] =?UTF-8?q?Sch=C3=A9ma=20GKS=20=C3=A9crit=20Erreur=20oscil?=
=?UTF-8?q?lation=20fin=20de=20l'onde=20de=20d=C3=A9tente=20D=C3=A9placeme?=
=?UTF-8?q?nt=20de=20l'onde=20de=20d=C3=A9tente,=20choc=20et=20discontinui?=
=?UTF-8?q?t=C3=A9=20de=20contact=20lrosque=20dt=20est=20modifi=C3=A9?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
---
src/language/modules/SchemeModule.cpp | 6 +-
src/scheme/GKS.cpp | 125 +++++++++++++++++++++-----
src/scheme/GKS.hpp | 1 +
3 files changed, 108 insertions(+), 24 deletions(-)
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 0be440945..76dfefd4a 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 09fa5f4ff..c3f296b73 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 cd3643003..e29ad2aaa 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
--
GitLab