Skip to content
Snippets Groups Projects
Commit 5eb27297 authored by clovis schoeck's avatar clovis schoeck
Browse files

Programme GKS à l'ordre 1 tourne et résout les équations d'Euler

parent d77b1239
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/language/modules/SchemeModule.cpp
+ 3
− 2
View file @ 5eb27297
......@@ -426,11 +426,12 @@ SchemeModule::SchemeModule()
[](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_U,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_E, const double gamma,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_E,
const std::shared_ptr<const DiscreteFunctionVariant>& tau, const double& delta,
const double& dt) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>> {
return gks(rho, rho_U, rho_E, gamma, dt);
return gks(rho, rho_U, rho_E, tau, delta, dt);
}
));
......
src/scheme/GKS.cpp
+ 120
− 48
View file @ 5eb27297
......@@ -17,15 +17,24 @@ class GKS
std::shared_ptr<const DiscreteFunctionVariant> rho_v,
std::shared_ptr<const DiscreteFunctionVariant> rho_U_v,
std::shared_ptr<const DiscreteFunctionVariant> rho_E_v,
const double gamma,
std::shared_ptr<const DiscreteFunctionVariant> tau,
const double delta,
double dt)
{
using Rd = TinyVector<MeshType::Dimension>;
const MeshType& mesh = *p_mesh;
const double pi = acos(-1);
const int delta = (3 - gamma) / (gamma - 1);
const double pi = std::acos(-1);
DiscreteFunctionP0<const double> tau_n = tau->get<DiscreteFunctionP0<const double>>();
CellValue<double> eta(mesh.connectivity());
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
if (tau_n[cell_id] == 0)
eta[cell_id] = 0;
else
eta[cell_id] = (tau_n[cell_id] / dt) * (1 - std::exp(-dt / tau_n[cell_id]));
}
DiscreteFunctionP0<const double> rho_n = rho_v->get<DiscreteFunctionP0<const double>>();
DiscreteFunctionP0<const Rd> rho_U_n = rho_U_v->get<DiscreteFunctionP0<const Rd>>();
......@@ -36,101 +45,163 @@ class GKS
DiscreteFunctionP0<double> rho_E = copy(rho_E_n);
CellValue<double> lambda{p_mesh->connectivity()};
// lambda.fill(0);
CellValue<Rd> U{p_mesh->connectivity()};
// U.fill(0);
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_flux_Euler(mesh.connectivity());
NodeValue<Rd> rho_U_flux_Euler(mesh.connectivity());
NodeValue<double> rho_E_flux_Euler(mesh.connectivity());
rho_flux_Euler.fill(0);
rho_U_flux_Euler.fill(TinyVector<1>(0));
rho_E_flux_Euler.fill(0);
NodeValue<double> rho_flux_Navier(mesh.connectivity());
NodeValue<Rd> rho_U_flux_Navier(mesh.connectivity());
NodeValue<double> rho_E_flux_Navier(mesh.connectivity());
rho_flux_Navier.fill(0);
rho_U_flux_Navier.fill(TinyVector<1>(0));
rho_E_flux_Navier.fill(0);
NodeValue<double> rho_node(mesh.connectivity());
NodeValue<Rd> rho_U_node(mesh.connectivity());
NodeValue<double> rho_E_node(mesh.connectivity());
rho_node.fill(0);
rho_U_node.fill(TinyVector<1>(0));
rho_E_node.fill(0);
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);
lambda[cell_id] = 0.5 * (1. + delta) * rho_n[cell_id] / (2 * rho_E_n[cell_id] - rho_U_2);
}
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
double U_2 = U[cell_id][0] * U[cell_id][0];
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;
rho_node[node_list[1]] = 0.5 * rho_cell_left;
rho_U_node[node_list[1]][0] = 0.5 * rho_U_cell_left[0];
rho_E_node[node_list[1]] = 0.5 * rho_E_cell_left;
}
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
for (CellId cell_id = 1; cell_id < mesh.numberOfCells(); ++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.;
double rho_cell_right = rho_n[cell_id] * (1. - std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0]));
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_node[node_list[0]] += 0.5 * rho_cell_right;
rho_U_node[node_list[0]][0] += 0.5 * rho_U_cell_right[0];
rho_E_node[node_list[0]] += 0.5 * rho_E_cell_right;
}
for (CellId cell_id = 1; cell_id < mesh.numberOfCells(); ++cell_id) {
auto node_list = cell_to_node_matrix[cell_id];
double rho_U_2_node = rho_U_node[node_list[0]][0] * rho_U_node[node_list[0]][0] / rho_node[node_list[0]];
rho_flux_Euler[node_list[0]] = rho_U_node[node_list[0]][0];
rho_U_flux_Euler[node_list[0]][0] =
delta * rho_U_2_node / (1. + delta) + 2 * rho_E_node[node_list[0]] / (1. + delta);
rho_E_flux_Euler[node_list[0]] =
rho_U_node[node_list[0]][0] / rho_node[node_list[0]] *
((3. + delta) * rho_E_node[node_list[0]] / (1. + delta) - rho_U_2_node / (1. + delta));
}
//##
//%%%%%%%%%%%%%%%
//##
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++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];
Rd F2_fn_left;
F2_fn_left[0] = (rho_U_2 + 0.5 * rho_n[cell_id] / lambda[cell_id]) *
(1. + std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0])) +
rho_U_n[cell_id][0] * std::exp(-lambda[cell_id] * U_2) / std::sqrt(pi * lambda[cell_id]);
double F3_fn_left = 0.5 * rho_U_n[cell_id][0] * (U_2 + 0.5 * (delta + 3) / lambda[cell_id]) *
(1. + std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0])) +
0.5 * rho_n[cell_id] * (U_2 + 0.5 * (delta + 2) / lambda[cell_id]) *
std::exp(-lambda[cell_id] * U_2) / std::sqrt(pi * lambda[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]]);
rho_U_flux_Navier[node_list[1]][0] += 0.5 * F2_fn_left[0];
rho_E_flux_Navier[node_list[1]] += 0.5 * F3_fn_left;
}
// 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(); ++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];
Rd F2_fn_right;
F2_fn_right[0] = (rho_U_2 + 0.5 * rho_n[cell_id] / lambda[cell_id]) *
(1. - std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0])) -
rho_U_n[cell_id][0] * std::exp(-lambda[cell_id] * U_2) / std::sqrt(pi * lambda[cell_id]);
double F3_fn_right = 0.5 * rho_U_n[cell_id][0] * (U_2 + 0.5 * (delta + 3) / lambda[cell_id]) *
(1. - std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0])) -
0.5 * rho_n[cell_id] * (U_2 + 0.5 * (delta + 2) / lambda[cell_id]) *
std::exp(-lambda[cell_id] * U_2) / std::sqrt(pi * lambda[cell_id]);
// double F3_fn_right = U[cell_id][0] * (rho_E_n[cell_id] + 0.5 * rho_n[cell_id] / lambda[cell_id]) *
// (1. - std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0])) -
// 0.5 * (rho_E_n[cell_id] + 0.25 * rho_n[cell_id] / lambda[cell_id]) *
// std::exp(-lambda[cell_id] * U_2) / std::sqrt(pi * lambda[cell_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 * 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];
rho_U_flux_Navier[node_list[0]][0] += 0.5 * F2_fn_right[0];
rho_E_flux_Navier[node_list[0]] += 0.5 * F3_fn_right;
}
// std::cout << "lambda " << lambda << std::endl;
// std::cout << "rho flux " << rho_flux_Euler << std::endl;
// std::cout << "rhoU flux " << rho_U_flux_Euler << std::endl;
// std::cout << "rhoE flux " << rho_E_flux_Euler << std::endl;
// std::exit(0);
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_Euler_sum = (rho_flux_Euler[node_list[1]] - rho_flux_Euler[node_list[0]]);
const Rd rho_U_flux_Euler_sum = (rho_U_flux_Euler[node_list[1]] - rho_U_flux_Euler[node_list[0]]);
const double rho_E_flux_Euler_sum = (rho_E_flux_Euler[node_list[1]] - rho_E_flux_Euler[node_list[0]]);
const Rd rho_U_flux_Navier_sum = (rho_U_flux_Navier[node_list[1]] - rho_U_flux_Navier[node_list[0]]);
const double rho_E_flux_Navier_sum = (rho_E_flux_Navier[node_list[1]] - rho_E_flux_Navier[node_list[0]]);
rho[cell_id] -= dt / Vj[cell_id] * (rho_flux_Euler_sum);
rho_U[cell_id][0] -=
dt / Vj[cell_id] *
(rho_U_flux_Euler_sum[0] + eta[cell_id] * (rho_U_flux_Navier_sum[0] - rho_U_flux_Euler_sum[0]));
rho_E[cell_id] -=
dt / Vj[cell_id] * (rho_E_flux_Euler_sum + 0 * eta[cell_id] * (rho_E_flux_Navier_sum - rho_E_flux_Euler_sum));
}
return std::make_tuple(std::make_shared<DiscreteFunctionVariant>(rho),
std::make_shared<DiscreteFunctionVariant>(rho_U),
std::make_shared<DiscreteFunctionVariant>(rho_E));
......@@ -145,7 +216,8 @@ std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, // rho
gks(std::shared_ptr<const DiscreteFunctionVariant> rho_v,
std::shared_ptr<const DiscreteFunctionVariant> rho_U_v,
std::shared_ptr<const DiscreteFunctionVariant> rho_E_v,
const double gamma,
std::shared_ptr<const DiscreteFunctionVariant> tau,
const double delta,
const double dt)
{
std::shared_ptr mesh_v = getCommonMesh({rho_v, rho_U_v, rho_E_v});
......@@ -165,7 +237,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, rho_U_v, rho_E_v, gamma, dt);
return gks.solve(p_mesh, rho_v, rho_U_v, rho_E_v, tau, delta, dt);
} else {
throw NormalError("dimension not treated");
......
src/scheme/GKS.hpp
+ 6
− 5
View file @ 5eb27297
......@@ -6,12 +6,13 @@
#include <tuple>
std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, // rho
std::shared_ptr<const DiscreteFunctionVariant>, // U
std::shared_ptr<const DiscreteFunctionVariant>> // E
std::shared_ptr<const DiscreteFunctionVariant>, // rhoU
std::shared_ptr<const DiscreteFunctionVariant>> // rhoE
gks(std::shared_ptr<const DiscreteFunctionVariant> rho,
std::shared_ptr<const DiscreteFunctionVariant> U,
std::shared_ptr<const DiscreteFunctionVariant> E,
const double gamma,
std::shared_ptr<const DiscreteFunctionVariant> rhoU,
std::shared_ptr<const DiscreteFunctionVariant> rhoE,
std::shared_ptr<const DiscreteFunctionVariant> tau,
const double delta,
const double dt);
#endif // GKS_HPP
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment